{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weighted Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "from statsmodels.sandbox.regression.predstd import wls_prediction_std\n",
    "from statsmodels.iolib.table import (SimpleTable, default_txt_fmt)\n",
    "np.random.seed(1024)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## WLS Estimation\n",
    "\n",
    "### Artificial data: Heteroscedasticity 2 groups \n",
    "\n",
    "Model assumptions:\n",
    "\n",
    " * Misspecification: true model is quadratic, estimate only linear\n",
    " * Independent noise/error term\n",
    " * Two groups for error variance, low and high variance groups"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "nsample = 50\n",
    "x = np.linspace(0, 20, nsample)\n",
    "X = np.column_stack((x, (x - 5)**2))\n",
    "X = sm.add_constant(X)\n",
    "beta = [5., 0.5, -0.01]\n",
    "sig = 0.5\n",
    "w = np.ones(nsample)\n",
    "w[nsample * 6//10:] = 3\n",
    "y_true = np.dot(X, beta)\n",
    "e = np.random.normal(size=nsample)\n",
    "y = y_true + sig * w * e \n",
    "X = X[:,[0,1]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### WLS knowing the true variance ratio of heteroscedasticity\n",
    "\n",
    "In this example, `w` is the standard deviation of the error.  `WLS` requires that the weights are proportional to the inverse of the error variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.927\n",
      "Model:                            WLS   Adj. R-squared:                  0.926\n",
      "Method:                 Least Squares   F-statistic:                     613.2\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):           5.44e-29\n",
      "Time:                        15:06:39   Log-Likelihood:                -51.136\n",
      "No. Observations:                  50   AIC:                             106.3\n",
      "Df Residuals:                      48   BIC:                             110.1\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2469      0.143     36.790      0.000       4.960       5.534\n",
      "x1             0.4466      0.018     24.764      0.000       0.410       0.483\n",
      "==============================================================================\n",
      "Omnibus:                        0.407   Durbin-Watson:                   2.317\n",
      "Prob(Omnibus):                  0.816   Jarque-Bera (JB):                0.103\n",
      "Skew:                          -0.104   Prob(JB):                        0.950\n",
      "Kurtosis:                       3.075   Cond. No.                         14.6\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "mod_wls = sm.WLS(y, X, weights=1./(w ** 2))\n",
    "res_wls = mod_wls.fit()\n",
    "print(res_wls.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OLS vs. WLS\n",
    "\n",
    "Estimate an OLS model for comparison: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5.24256099 0.43486879]\n",
      "[5.24685499 0.44658241]\n"
     ]
    }
   ],
   "source": [
    "res_ols = sm.OLS(y, X).fit()\n",
    "print(res_ols.params)\n",
    "print(res_wls.params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the WLS standard errors to  heteroscedasticity corrected OLS standard errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=====================\n",
      "          x1   const \n",
      "---------------------\n",
      "WLS     0.1426  0.018\n",
      "OLS     0.2707 0.0233\n",
      "OLS_HC0  0.194 0.0281\n",
      "OLS_HC1  0.198 0.0287\n",
      "OLS_HC3 0.2003  0.029\n",
      "OLS_HC3  0.207   0.03\n",
      "---------------------\n"
     ]
    }
   ],
   "source": [
    "se = np.vstack([[res_wls.bse], [res_ols.bse], [res_ols.HC0_se], \n",
    "                [res_ols.HC1_se], [res_ols.HC2_se], [res_ols.HC3_se]])\n",
    "se = np.round(se,4)\n",
    "colnames = ['x1', 'const']\n",
    "rownames = ['WLS', 'OLS', 'OLS_HC0', 'OLS_HC1', 'OLS_HC3', 'OLS_HC3']\n",
    "tabl = SimpleTable(se, colnames, rownames, txt_fmt=default_txt_fmt)\n",
    "print(tabl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate OLS prediction interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "covb = res_ols.cov_params()\n",
    "prediction_var = res_ols.mse_resid + (X * np.dot(covb,X.T).T).sum(1)\n",
    "prediction_std = np.sqrt(prediction_var)\n",
    "tppf = stats.t.ppf(0.975, res_ols.df_resid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "prstd_ols, iv_l_ols, iv_u_ols = wls_prediction_std(res_ols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot to compare predicted values in WLS and OLS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAFnCAYAAACRo/HLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd8zdf/wPHXJzsRhNiE2FGjRtBrVFSNVrVGq61dVJXao9Whg19J7VJVm1KtWh2K9ktjxohN7RkxI4LIvPee3x9H7BHJTe5NvJ+Px32oez+59yTR+77nnPd5vw2lFEIIIYTIeE72HoAQQgjxtJIgLIQQQtiJBGEhhBDCTiQICyGEEHYiQVgIIYSwEwnCQgghhJ1IEBZCCCHsRIKwEEIIYScpCsKGYeQ3DGP9zf8uYRjGasMwdhmG8Vn6Dk8IIYTIuh4bhA3DyAXMAbLdvOsDYKhSqjLQ2DCMvOk4PiGEECLLcknBNRbgTeC3m3+/DFQyDOMo4A5EP+qL8+TJo/z9/dMyRiGEECLT2L59e6RSKkUT1McGYaXUNQDDMJLvWgn0BooAawDzvV9jGEY3oBtA0aJFCQsLS9HAhRBCiMzOMIxTKb02NYlZHwGdlFKfAJ5Aw3svUEpNVUoFKqUC8+aV1WohhBDiQVIThIsDfoZheABVAWnDJIQQQqRCaoLw50AIcAkIRy9JCyGEEOIJpSQxCwClVNDNP5cDy9PyoklJSZw5c4b4+Pi0PI1D8vDwoEiRIri6utp7KEIIIRxcioOwLZ05c4bs2bPj7+9/Z8JXpqeU4vLly5w5c4bixYvbezhCCCEcnF0qZsXHx+Pr65ulAjDoDHJfX98sOcMXQghhe3YrW5nVAnCyrPp9CSGEsL1MUTt62c4Iao9cQ/GPllN75BqW7YxI0/N16tSJypUrExgYyLRp0x56XUhICCdPnkzTawkhhBAPY5c94SexbGcEQ5bsJS7JAkBEdBxDluwFoHmVwql+3kmTJlGuXDmeffZZatasSaVKle67JiQkhKCgIKTilxBCiPTg8DPhUasO3QrAyeKSLIxadSjNz+3r60vTpk1ZtmwZderUoW7dunzyyScAvPPOO8yePZu+ffvStm1bAM6ePXvfdUIIIRyTrVdR04PDB+Gz0XFPdP+TSk4QGzlyJCtWrOCPP/4AYNasWXTq1Inx48czf/58ACIiIu67TgghhONJXkWNiI5DcXsV1dECscMvRxfy8STiAQG3kI+nTZ4/KioKPz8/Ro4cibe3N9evX3/otS4uLnz55ZePvU4IIYR9PWoVNS1bmbbm8DPhQY3L4unqfNd9nq7ODGpcNs3PHR0dzYoVK9i1axdDhgxh+vTpd2U3e3p6EhsbC+gzwGPHjn3gdUIIIRxLeq+i2orDz4STP7GMWnWIs9FxFPLxZFDjsmn+JNOrVy/c3d0JDg7GMAy6d+9O3rx58fLyIiIigsKFC9OqVSs6d+7MV199xbx583jllVceeJ0QQgjHkt6rqLZiKJW+/RcCAwPVva0MDxw4QLly5dL1de0pq39/Qgjh6O49WQN6FXVEy4rpvhxtGMZ2pVRgSq51+JmwEEII8aTSaxXV1iQICyGEyJKaVynscEH3Xg6fmCWEEEJkVRKEhRBCCDuRICyEEELYyVMbhCdNmkRQUBCenp4EBQWxdOlSew9JCCHEU+apDcIffPABISEhFC5cmJCQEFq0aGHvIQkhhMhgRy4fIb2P6j6K3bOj+/aFXbts+5yVK8P48U/+dUFBQVSvXp09e/awatUqvvjiC4KCgggKCmL27NkAtG7dmg4dOnDx4kUqVqzId999Z9vBCyGESHcno0/Se0Vv/jj8B2s6rKF+8fp2GcdTOxN+kM2bN2MymVi1atVDr5k6dSoVKlRg3bp1nDt3jj179mTgCIUQQqSWUooLMRcAyOmek70X9/Jl0JdUyn9/K9uMYveZcGpmrOmlQoUKtGzZ8oGPxcXF4enpyaFDh9i0aRMhISFER0cTERHxwF7EQgghHIPZauaXfb8QvDEYFycXtnfbTi7PXBztdRRnJ+fHP0E6kpnwHby9ve/6u5ubG5cuXQJg5cqVAJQtW5a+ffsSEhLC8OHDKVq0aIaPUwghxOPFJsUyaeskSk8sTbul7bAqK32f64tC7wHbOwCDA8yEHdmrr75Kjx49WL16Nb6+vgC8++67vPPOO8yaNYscOXLw008/2XmUQgghHuSXfb/Qa0UvTEVMfNvkW5qWaYqT4VhzT2ngkA6y+vcnhBCO6My1M4wLHUfZPGXpVq0bCeYEtp3dRp2idTJ0HNLAQQghxFPjwKUDjNo0inl75mFVVvqb+gPg7uKe4QH4SUkQFkIIkWl9/u/nDFs3DA8XD96r9h4Dag3A38ff3sNKMQnCQgghMg2lFH8f+5vKBSqT3zs/dYrW4dPnP6VXjV7kzZbX3sN7Yo61Qy2EEEI8gNlq5ud9P1N1alWazG/C1O1TAWhYsiFf1f8qUwZgkJmwEEIIBzd1+1SCNwZz/MpxAvIEMOu1WbSp2Mbew7KJpzIIx8TE0L59ey5dukTJkiXx8/MjICCAdu3a3brmxo0btGvXjqioKIoWLcrcuXMxDMOOoxZCiKdHbFIsXq5eAPxz/B/yeuVlTKMxvFr2VYc7ZpQWWec7eQITJ06kdOnSbNiwgYSEBBYuXHjfNT/++CMmk4m1a9fi7u7OvceshBBC2N7Z62cZ9PcgCo0pxIFLBwCY/dpsQruE0jygeZYKwOAoM+GgoPvva90aevSA2Fh4+eX7H+/USd8iI+H11+9+LCTkkS+3ZcsWunbtCkCdOnUoWLDgfdcULlyYOXPm0KJFC6ZPn56ib0MIIUTqHIo8xKhNo/hxz4+YrWbeLP8mbs5uAGRzy2bn0aUfxwjCGez69etky6Z/qV5eXly7du2+a5o1a0ZcXBwtW7akfv36jBs3Dmdn+5c4E0KIrCYmMYZqU6thURa6VunKgFoDKJGrhL2HlSEcIwg/aubq5fXox/PkeezM9145cuQgJiYG0Hu/OXLkuO+aI0eO0KRJE1q1akW7du2YN28eHTt2fKLXEUIIcT+lFP8c/4c/D//JhCYT8HbzZkGrBdQoXIP83vntPbwMlaLFdcMw8huGsf6e+/4wDKNy+gwrfdWsWZOQm4F7/fr1/PXXX/ddM336dJYuXYqzszMVKlQgPj4+g0cphBBZi8Vq4Zd9v1BtajUaz2vM4gOLOR9zHoBmZZs9dQEYUlA72jCMXMACIJ9SqurN+9oC1ZVSfR/3Ao5YOzo5O/rChQuULl0aPz8/fvrpJ3Lnzg1Ap06daNmyJW3btkUpRc6cOVmwYAFeXl4pen57f39CCOFo9l/cz2s/v8axK8co61uWQbUG0a5SO9xd3O09NJuzde1oC/Am8NvNJ88NjAG+NwyjvlLq31SP1E68vb1ZunTpXfcNHz78vuv+/TfTfWtCCOEwouOjOXHlBFUKVqFErhKU8S3DqIajeC3gtSyX5Zxajw3CSqlrwJ1nZPsBvwI/ACMMw8iulPr9zq8xDKMb0A2QfrtCCPGUOXf9HOM3j+f7sO/Jmy0vR3odwdPVk7/a3r/1ZzfJq8B2rv+Qmo8iVYDvlFLngYVA0L0XKKWmKqUClVKBefNmzlJiQgghnsyxqGN0+6Mb/hP8GR06mqZlmrLojUWONeu1WGDRIqhZEzZutPdoUpUdfRQoARwEAoFTNh2REEKITMWqrDgZTuy+sJu5u+fSpUoXBpgGUDJ3SXsP7ba4OJgzB0aPhmPHoGRJuHlKxp5SE4S/AaYbhvEJEAu0tO2QhBBCOLrkY0bBG4OpW7QuXwR9QfOA5pzqe8rxspyVgurVYf9+/WdwMDRvDg5Q+yHFQVgpFXTzz7PAA0pYCSGEyOrMVjOL/1tM8MZgdp7fSUHvgrzxzBsAOBlOjhOAT56EWbNg6FAdbD/9FAoWhOeft/s+8J0caKE+45QtW5aoqCh8fX05e/Ys9evXZ968eXddc+PGDVq0aEG9evVo3749jzvKJYQQT4Oey3vy1uK3iE2KZXqz6Zzoc4Lugd3tPazbdu6ENm2gVCn4+mvYsUPf/9ZbUK+eQwVgeEqDcJ48edi+fTtXr15l//79FCtW7L5rpIGDEELoY0Yj1o/gaNRRALoHdmdJ6yX81/M/ulTt4jjnfC9dgkaNoGpV+PNP6NsXjh/Xy88OzCHKVgbNDrrvvtblW9Ojeg9ik2J5ef79q9+dKneiU+VORMZG8vrCuxs4hHQKeeTrFStWjLVr11KnTh3Wrl37wCAsDRyEEE+zO48ZXU+8jperF32e60OVglWoUrCKvYenmc1w6BCULw+5c0NCAowcCe+9Bz4+9h5dijhEEM5o/v7+rFu3jiZNmrBy5UpOnTrF//3f/911jTRwEEI8rXr91YupO6ZitpppXb41g2sNdpzACzqrecYMGDcObtyAU6d0n4G1a+09sifmEEH4UTNXL1evRz6exyvPY2e+9ypWrBijR49mzJgxDB06lFq1at13jTRwEEI8TfZf3E/5fOUBnWDlkMeMLl6Eb7+FyZPhyhWoUwcGDQIPD3uPLNWeyj1hf39/8uXLR/ny5bFYLPj7+993jTRwEEJkdUopVh1dxQtzXqDC9xUIDQ8FYMJLE5jcdLLjBGCrVf95+LBOtgoKgk2bYP16ePVVcMq8ocwhZsIZrVixYpQsWRIvLy+KFCmCn58fQ4cOZfz48YBu4NCnTx/atm3LrFmzbjVwEEKIrCD5mNHIjSPZdX4XhbIXYnTD0VTIV8HeQ7tbaCiMGgWFC8PEiVC7tk62esDEKbN6KoNwQEAAa2/uHYSHhwPctycM0sBBCJE1xZvjeX/5++TLlo8Zr86gbcW2jpPlbLXq7OZRo2DDBp1gNWCAfswwslQAhqc0CAshxNPkStwVJm+bzN/H/+bfjv/i7ebNpi6bKONbxrHqOgN89pleci5aFMaPhy5dwNvb3qNKNxKEhRAii4q4FsG4zeP4YfsPxCTG8FKpl7gSdwVfL18C8gTYe3hadDRMmQIvvgiBgdCpkz5y9MYb4Opq79GlO7sFYaXUne0RswyprCWEcARbzmyh7qy6WJSFtyq8xeBag3m2wLP2HtZt4eF6pjt1qj5ylJSkg3Dp0vr2lLBLEPbw8ODy5cv4+vpmqUCslOLy5ct4ZOJ0eSFE5rXlzBYirkfQslxLqhWqxuDag+lSpQvFcxW399Du1qePPmakFLz5pj5mVLmyvUdlF3YJwkWKFOHMmTNcunTJHi+frjw8PChSpIi9hyGEeEoopVh5dCXBG4NZe2otAXkCaBHQAhcnF4a/MNzew9OU0seJ6tTRx4ny54eePaFfP3hAxcKniV2CsKurK8WLO9gnMyGEyGT+PfEv/Vb1Y/eF3RTJUYSxjcbybrV3HWeF0WyGRYt0pvOOHfD779CsGXz8sb1H5jAkMUsIITKR2KRY4s3x5PbMjWEYJFoSmfXaLNpUbIObs5u9h6clJsIPP8DYsbqlYJkyMG0aNGxo75E5HCO9E4kCAwOVdCASQoi0uRx7me+2fcfErRN5u8LbfPvStyilUCjHOWaUlKQzms1m3UqwcGG935vJq1o9KcMwtiulAlNyrcyEhRDCgYVfDWds6Fim7ZjGjaQbNC3dlDfLvwmAYRgYOMDS85EjMGYMrFoFBw7oWs7btkHevPYemcOTICyEEA7s85DPmbt7Lm0qtmFQrUFUzF/R3kO6bfNm+OYbWLYM3NygQwfd1cjDQwJwCslytBBCOJCNpzcSvDGYz+t9TrVC1Th99TRKKYr5OFgW8bZtUKMG5MoFPXpAr14661nIcrQQQmQmVmVl+eHlBG8MZmP4RnJ75uZE9AmqFapG0ZxF7T08LT4e5s2Dq1d1LefAQJg7F1q0yNJlJdObBGEhhLAjpRTPz3qejeEbKZqzKBOaTKBLlS5kc8tm76FpV67A99/rPr4XLkDdutC/v26m0L69vUeX6T096WpCCOEgYhJjmLVz1q3yvW0qtuHHFj9ytNdRetfs7TgBeMEC8PODTz7RFa3+9z9Yu1YHYGETMhMWQogMcunGJSZunch3274jKi6K0r6lqVO0Dj2q97D30G7buROyZdNneytUgJYtYeBAqFTJ3iPLkmQmLIQQ6exq/FV6/dWLYuOLMWzdMJ4v9jybOm+iTtE69h6aphSsXKk7GVWtCsn91StW1Pu+EoDTjcyEhRAinUTHR+Pj4YOXqxcrjq7grQpvMajWIMrlLWfvod32668wbBjs3QuFCkFwMLz3nr1H9dSQICyEEDaklGLdqXWM3DiSvRf2cqz3Mdxd3Pmv53+OU1by+nWd0WwY+qiR1QqzZkGbNvq8r8gwshwthBA2YLFaWHJgCc/NeI6gOUFsP7ud9wPfx6IsAI4RgCMi4MMPoUgRXd0K4MsvYc8e6NRJArAdyExYCCFsYO2ptbRa2IoSuUow+eXJdKrcCU9XT3sPS9u3D0aPhp9+AosFXn9dZz0DeDrIGJ9SEoSFECIVriVcY0rYFJRSfFjnQ+r71+fPt/+kSakmODs523t4t1ks8PLLcPkydO+ue/hKK1mHIUFYCCGewPmY80zYPIHJYZO5lnCNFgEtbp33bVqmqb2HpzsYLV6sq1stXqyXmBcuhNKlwdc3w4axbGcEo1Yd4mx0HIV8PBnUuCzNqxTOsNfPLCQICyFECs3cOZP3l7+P2WqmVblWDK49mMBCKSoRnP5iYmDmTBg37nYP31OndPB97rkMHcqynREMWbKXuCS9Hx4RHceQJXsBJBDfQxKzhBDiEbZFbONo1FEAqhWsxjuV3+HQB4dY+MZCxwnAx45B0aLQp4/u4btsmW4pWLq0XYYzatWhWwE4WVyShVGrDtllPI8SGh7KiPUjCA0Ptcvry0xYCCHuoZTin+P/ELwxmDUn1tClShemvzqdZws8y5RXpth7eNrBg7B/P7RqBSVKQMeO0Lo1mEz2Hhlno+Oe6H57CQ0Ppf6c+pitZtyc3VjdYTUmv4z9+UkQFkKIOyw9sJRh64ax8/xOCnoX5JsXv+G9QAcpXqEUbNigM51//123Dnz1VXB11cvQDqKQjycRDwi4hXxsm4md1n3nYeuGkWBJACDRkkjIyZAMD8IpWo42DCO/YRjr77mvgmEY/6TPsIQQIuPEJcWR3Ft93al1xCbFMr3ZdE70OcGg2oPI4Z7DziMENm/We7vPPw8bN8LQofp8r6urvUd2n0GNy+LpeneGuKerM4Mal7XZayTvO0dEx6G4ve+8bGfEQ78mwZzAjB0z2HJmCwBdq3bF1ckVZ8MZN2c3gvyDbDa+lHrsTNgwjFzAHCDbHfcZwFjA8X77QgiRQlFxUXy39Tu+3fotP7f6mQYlGjD8heGMaTwGJ8MBUmZiY3XCVb58Oss5KgomT9ZLz15e9h7dQyXPRtMzO/pR+853vk5oeCgrj64kMjaSpQeXci7mHL1q9KJmkZq0LNeStZ3WEnIyhCD/oAyfBUPKlqMtwJvAb3fc9w7wL9A4PQYlhBDpKfxqOGNDxzJtxzRuJN3g5dIvk9szN4BjtBG8eBEmTdIBt2lTmDNHN1Y4dAicHODDQQo0r1I4XTOhU7LvHBoeSr3Z9UiyJgEQWCiQOc3n8GKJF29dY/Iz2SX4JntsEFZKXQMwbvaPNAzDF2iHDsAPDMKGYXQDugEULVrURkMVQoi0syordWbVIeJaBG9XfJvBtQZTMX9Few9LO3wYxozRQTchQe/3dut2+/FMEoAzwqP2nQ9cOkBp39KEnAzBbDUD4Gw40zKgJQ1LNszooT5San6jI4EhSqmkh12glJqqlApUSgXmzZs39aMTQggb2HB6A11+60KSJQknw4mZr87kWO9j/NjiR8cIwDf3o5k8WQfgjh31EaPffoPate07Ngd1776zQqFcD5DoM4JnJj/DkgNLCPIPwsPF4/49X7MZfvkFIiPtM/g7GMnJCI+90DBClFJBhmEcBs7evLsyMEkp9enDvi4wMFCFhYWlfaRCCPEErMrKn4f/JHhjMJvCN+Hr6UtIpxAq5Ktg76FpFosOsqNHw/Dh8MILehlaKZ31nI6ySjWrZTsj+OyvxZyIX4pyOUksJ/H19KVXjV70rNGTPF55CA0Pvb3nm7M8zJgBEyboQibjx+uz1TZmGMZ2pVSKDpE/8RElpVSZO14o5FEBWAgh7CHiWgQNf2zIgcgDFMtZjIkvTaRzlc54uTpAMlNsrJ7tjh0LR4/qOs6xsfqxfPnS/eUzSzWrx31QUEqRP89pjvERcS5xGBj0f64/w14Ydtfv2eRnwlS4Jnz0EfzQBK5dg7p1dSBu1swe39pdUhyElVJBKblPCCHs4XrCdXZf2E2donUomL0gFfJV4JO6n9C6fGtcnR3kIIdS+pjR3r1Qvbqu6dyyJThnXMOHlGYV29OjPigElcvGlLApzNszjzeeeYNESyIAToYTebzy3P1B6+RJ8PfXe+mHDsFLL8GAAfpn7yCkWIcQIlO7EHOBb7d8y+SwyQBE9I/Ay9WLhW8stPPIbjpyRC+BDh8OLi7w+eeQN6+ejd1MeM1ImaGa1YM+KFxPukDP5TO5tnIFMYkxNCrZiIr5K+Lm7EaiJfH2nq/VCn/9pRPc1q27vdqwZEmGfthJKQnCQohM6fTV04xYP4JZu2aRaEmkZbmWDK492DGWnAE2bdL7vcuW6TO+rVrpGVirVnYdVkZVs0qL5A8ECU4HiHfai4u1MJFu34BZ0eaZtxhUaxCVC1QGYHWH1XrPt5AJ04q9MK6zLulZpAgEB9/uHOWAARgkCAshMpkkSxKuzq6cjznPzF0z6fhsRwbWGkgZ3zKP/+KMcPEiNG8OoaGQOzd88gl88EG6J1ul1KDGZe9a6gXbV7NKq4I5PTgYs5grrtMAKwYuZDe/TOlsbzG/ZZu7rjUVeU6f8z17Fj5oBJUqwU8/weuvO2Q1sXtJEBZCODylFKtPrCZ4YzBFcxRlxmszqFG4BhH9I8jjlcfew9OJVfv2QY0akCcP5MoFEyfCO+9ANgco/nGHjKhmlVoWq4XfD/3O5WzDuZK4AxRggFJmPJx8+axJvdsXHzigk9vOnIEVK6BQIb3XXqaMXZb5U0uCsBDCYZmtZhb/t5hvNn3DjnM7KOhdkKalm9563O4B+NIl+O47fbNYdEDw8oLly+07rsdI72pWqbH+1Hq6/tGVw5cPUyJXCV4r0Z0/js/Eqsw4Ga70r9eC5pULwZo1er/3r7/Aw0OfqU5K0rPeso4zm08pCcJCCIf1RcgX/N/6/6Osb1mmN5tOu0rtcHdxt/ewdNbtyJH6qFF8vK5sNXAgeDrOvmpmsOroKlafWE2LgBYU8C6Aj4cPC19fSMtyLXF2ciY0vMPddZ1nz9arC3nzwpdfwvvv6//OxCQICyEcxpW4K3y37TsaFG+Ayc9E16pdqVawGq8FvGb/hgpK6VKSHh5w4YIOCB06QP/+EBBg37FlMuFXwxn8z2B+3v8zAJO2TmJ1h9Vs6brlrutM2cthCl0PFyLAD53UZjZDu3b695AFSBAWQjgEpRR1ZtXhv0v/YVVWTH4m/H388ffxt+/ALBad4Tx6NFSoANOmQc2aOhEod277ji2T2XdxH6M2jeKnvT9hsd5ODLuvl++JE7qYxowZuotU9+460Sp7duja1U6jTx9SDVwI4RDOXDvDf5f+Y1TDUQytN9Tew4EbN/Reb5kyOgBcvKgTr5JJAH6s0PBQvl7/NZtObwJg4f6FLP5vMR9U/4BFrRfh6eJ5f13nzz6DUqX0z755c9ixA77/3n7fRDqTmbAQwiGEndU15usWrWvnkdz02Wcwbpye9X7zjQ4IDnrW1BFtOL2BBnMb3CqkEdIxhP6m/vR9ru+ttpGrO6wm5MQagk47Y8peTn9h5cp6f71XL33WN4uTICyEcAhhZ8NwcXLh2QLP2mcABw/qIy8dO+rORX366D3IWrUy1ZEXe4tLimPu7rl8vObjWyUlzVbz3cvNADduYPp9B6Zxs+DYMZjorc9Tt2pl94ImGUmCsBDCIbQu35rSvqXxcMnAhBuldGnDMWPgjz90sk+1ajoIFyumb+KJNJjbgNAzoQTkCSAmMQaL1XL3crPVCkOH6raNV67oWtojR+qVhqeQBGEhhEN4tsCzGT8LfuUVfd40Tx5d07lHjwzpZJTeMqJVYXKLwIA8AWyN2MrQekPxdPXkk7qfkM0tG/WK1WPzmc23jxgZfvoLnZwgLAzq19fNFGrVsum4MpsU9xNOLeknLIR4nEs3LhF6JpQg/yByuOdIvxeKiYGff9ZnTZ2d4Ycf9Gy4QwddZCMLuLcDEeiylCNaVrRZIA4ND6X+nPokWBIAcDacWdVuFQ1KNLj7QqXg77/1SsO//8Lx4+Dnp48ZuWTdOeCT9BOW7GghhN2FnAzhtZ9f42jU0fR5gXPn4OOPdQB4910ICdH3v/eePv6SRQIwPLpVoS1cT7hOx2UdbwVgA4MBpgF3B+CEBJg1S9dxbtJEl/T86it9xAiydAB+UhKEhRB2t+3sNtyc3aiQr4Jtn/jqVejcWe/tBgfDiy/qxgoNGjz+azOp9GhVaLFa2HV+FwDebt7ky5YPVydXnA1nPFw8aB5wcz83eWX14kXo1k0ntM2erc/9DhkCPj6pHkNWJR9HhBB2F3Y2jGfzP4ubs1van0wpOH1aB15vb9i+XQeEfv2gZMm0P7+Ds1WrwtDwUP45/g/XE66z9OBSzl4/y+l+p8njlYcNnTfc2hMO8g/ClJAXevbUBUyWLtUrDjt26OImkln+SBKEhRB2ZVVWtp/bTtuKbdP2RElJsHCh3n8MD4dTp/Qy844dT9X5Xlu0Klx1dBWvLHgFs9UMQLk85ZjXch65PHLdusZU5DlMpyzQ+xv47TfdQKFdu9v7vRUr2u6bysJkOVoIYVdHo45yLeEa1QtVT90TXLumA2+JEjoIxMfrIy/J+45PUQAG3SFpRMuKFPbxxAAK+3imOCnLqqwArDy68lYAdjKcaFep3a2mCqCTv/7vjQ+hbl2urlrNoS699YeeGTNkv/cJyU9LCGFXpXKX4mDPg+TN9oTdcJTSS527d+sKS/VgJ9NLAAAgAElEQVTr62znJk30MZin2JO2Ktx9fjejNo1CoZjfcj6ty7dmyvYpJFmScHN2o75/fbh+HWbMYDM5GXK5IK5FqhHXqAeLy78A2bIx4pyF5gXS8ZvKoiQICyHsyslwomyeJ+gDu3Onnvn6+uoi/3Xq6GbuFWyc1JXFbTq9iZm7ZrLv4j62RGzB282b9wPfRymFyc/Emg5r9J6v1zOYJi6FqS/B1aucr9GUuPrvE+fhzbwqL+snu5l97Wg9ih8lMRH27IGtW3VTrC+/tM84JAgLIezqm43fUC5POZqVbfbwi5SCVat0J6PVq3XCVe/e+jHDkAD8hELDQwmaE0SSNQmA7tW683WDr8nleceer58J07QVMOJ1/fN//XXo359+Sy498DnTkn2d3qxWOHpUB9zk286dOhADFC6sS4XbYyVdgrAQwm4sVgtfrv2SLlW6PDoIf/45DBum3y2Dg3W281N63CW11bBik2KZvWs2z+R9htDw0Fv7v86GM0VzFtUB2GrVFcSefx5y5IBy5XQjhd69wd8fgEJr1tgk+zo9Xb4MW7bA5s36tm0bREfrx7y8IDBQf0s1auhb0aL2S+KWICyEsJuDkQeJTYq9Pynr8mWYMgUaN9bvmO3b65aCrVuDmw2OMWVS91bDioiOY8iSvQAPDcSXYy/z3bbvmLh1IpGxkfQI7EG7Su1wc3a71eEoqKAJpk7VXaMOHtRtBHv0gLff1rc72CL72pbMZr0bkRxwQ0PhyBH9mJOTTtJ+4w3dDKtGDf25InnGu2xnBG0WpG95z8eRICyEsJvk9oWBhW5W+Dt2TAeCWbMgNlZPTwIDoXRpfXvKPaoa1r3BIzQ8lC9CvmDtqbUkWBJoWropg2sPpm7RuhiGcbuN4JrjmExvQGQkVK0K8+frqPUQya+T3rWpHyYqSgfbjRth0ya9tBwbqx/Llw9MJl2f5bnn9D8db+8HP09qPtCkBwnCQgi7CTsbhrebN2V8y+hykslHXNq1g/79Za/3HimphrXr/C6uxl/lpfkvEW+Ox8lw4scWP9KuUrvbX3Dhgt7z9TPBFy/oyNW/P9Srl6J12SfNvk4tpeDwYR1sk4PugQP6MWdnqFIFunTRPSCee07XZ0npsvKTfKBJTxKEhRD2YbFw7uhOqhWsps+flioFH32ke8oWKmTv0Tmkh1XDKpjTg9XHVxO8MZh/jv9Dq3KtSLQkotBlJMOvhuuIFhKiM8v/9z9dSrJgQVi50mGW+JOSdMLU+vWwYYO+RUbqx3Ll0sG2XTv9Z/XqkC1b6l8rPcp7poYEYSFExrpxQ9cTHjeORceOkfjnb/r+Dz+067Ayg+T92GjLPuKd9uJurYCzSxQXPJfz4o97KeBdgJENRlKlYBX+OvLX7T3fQwnQp5qOcHnz6mYWHjf7NtsxAMfE6KXl5KC7efPtpeWSJXWnydq19a1sWdse/7ZVec+0kiAshMgYN27oSlaTJ+uNvZo1YeRI3Jo0tffIMo3mVQpzMGo7n2z4FKtKwslwJZ9XQVxc3JjWbBrtKrXDw0UH19UdVutzvu5lMJla6yg2bZqeSiYH4AwWHa2XldeuhXXrdFlvs1kH18qVoWtXfey7Th09SU9PjpJgJkFYCJG+rl6FnDnB3V0n/Tz/PAwcyC85w/lp34/MTWpITuec9h5lphAZG8k/EROwkgAGGIaZTlXe4v8a/B9Oxs1p4smTMGECpvPnMS1YoO/btk1HuQyuJBYZqWe5yUF31y69Ku7mpjOVP9SVLzGZ9ImojGTvBLNkEoSFELanlG7iPnq0fuc9cUIH4b17b23k/ftnd9aeXEt29+x2HqzjO37lOGNDxzJz50zizHE4GU4YGLg5u/Fq2Vd1AN6yRe/3Ll6sg+3bb99uplC1aoaMMypKB9yQEP3r36uTjfH01IH288917lfNmvo+e8uoBLNHkSAshLCd5E5GycE3Xz6daGU26yB8RyZN2NkwAgsF3p7BiVvubBNYJEcRykwsc6uRwsBaA7kaf/V2G0E/E8ycqdOEc+bUdbR79YIiRdJ9nNHReoabHHR379afvzw99T7uW2/poFu9usPkfjkcCcJCCNvZuFHvOZYr98j9xwRzAnsu7KG/qb8dBvlwqa1GZUubTm/ihbkvkGhJxMPFg9UdVjO56WSalm5K4Rw3x3LjBqbfd0BMNPgBr72mGyx07gzZ029lIS5OJ1CtWaOrh27frotsubvrjOUvv9R9NKpX1/eJx0tREDYMIz+wSClV1zCMosBcwAocBd5TSql0HKMQwlGdOqWbKHh5wfDhetqzZo3+8xH7j3sv7iXJmnS7SIcDsHfxBrPVzK/7f2XA3wNIsCQAkGhJJORkCEPqDtEXnTsHkybpamJRUfps9Usv6WYWffrYfkxmvZ28erW+bdqk6y27uOgl5U8/1UH3uefsluuV6T02CBuGkQuYAySvI70HvK+UOmAYxgqgIrAn/YYohHA4YWF6//HXX/Xfu3TRfxqGfld+jARzArX8aqW+h3A6sGfxhvWn1tNhWQdORp+kWM5iuDq5YlVWfbzIP0hf9MUXMGKEXvJv0UIX16hd26bjUAoOHYJ//tG3kBA9wQZ49lm9s9CggU6mSscJ91MlJTNhC/Am8BuAUuqTOx7zBSLTYVxCCEc1YoQ+Z5o9O/TtqyvhFy36RE9Ru2htNnbemE4DTJ2MLN4QGh7K8iPLqZSvEq0rtMbfx59iOYsxockEXinzClvObCHk5L8EXfTClLuS/qLixXXjir599SFaG7l0Sc9y//5b1/AID9f3lyihc7tefBGCgvTxYltyhKV/R/DYIKyUugZg3FMLzDCMN4H9Sqmz6TM0IYRDiIuDuXN1emulSvDqq3rDr2vXVJ8rsVgtukqWA8mo4g2L/lvE24vexqzMOBlO+OX0w+RnIqRTiL4gIQHT/w5iGrsA9u2DKZ7w3nvQsaO+pVFiot66X7VKB96dO/X9Pj56lvvJJ9CwoQ7C6cXeS/+OJFVpiYZhlAAGAn0f8ng3wzDCDMMIu3Tpwb0nhRAO7uJFfaakaFHo3h1+/lnfX768XgpNZQCOTYol58icfL/texsONu0GNS6Lp+vdHwxsWbxh1/ldvLXoLVr/2hqzMgNgYBByMkRfYLHA11/rloGdO+ul/dmzoVOnNL2uUrqr0KRJ0KwZ5M4NL7wAY8fqX+Hw4fp0U2QkLFqk4316BmB49NL/0+aJs6Nv7hEvADorpa4+6Bql1FRgKkBgYKAkbQmR2QwaBBMnQkKCfuceMEAX2bCB3ed3cyPpBoWyO1Z96PQo3qCUQqFwMpz4+9jfrDi6gjYV27D4wGKSLEl6zzdXFX2xs7Pu5fvss3rl4cUXU93kNiZGLzGvWqVLQ584oe8vWVJPpps00UvM9trXdYi6zWaz/pnbq5HwTak5ovQRUBSYeHOJ+nOl1FqbjkoIkbGU0o1YTSb9puTpqWdg/frpcoc2dF/7Qgdiq+IN60+vZ9LWSew4u4Mvgr6gbaW29Kjeg/eqvUdOj5z0PN2DkHVzCFpxENPXrXSUzJdPZ0OlooqFUrq70IoVOo6vX6/zt7Jl00vMAwfq1sw23EpOE7vWbY6OhunT9YfMmTP1D8iOUhyElVJBN//8EJBK60JkBUlJOsN5zBjYsUNPmxo3hq++SreXDDsXRgHvAg43E7aFG4k3+HTNp0zYMgGFwsAg4noEAN5u3nr2tXAhpjFjMG3dqteG+/W73WX+CQJwTIw+DbZihb6dOqXvr1BB52699JJOnnbEIhl2qdt8+DB8+61e4r9xQx+jc4CyXVKsQ4inUVwcfPedPuN75gwEBMDUqTZbcn6U5EpZ9yZ7ZoT0zsh9af5LrD+9/tbfnQwnLNY79j5PnNBlpEqW1D//jh2fqB/fsWOwfLm+hYToJCtvb71y/fHHOvD6+dns20k3GV632WLRR+ciI3XKd58+uhmxAzDSu85GYGCgCgsLS9fXEEKkUHy8rqqQmKizb8qU0fu9L72UIcX9lVIEbwymRK4StC7fOt1f7073ZuSCnn2NaFkx1W/+x6KOMXHrRL6q/xU53HOw5sQajkYdpe/KvrqNoJMrq2NaYrrgBrNm6S8KDdXdC5wfnx2elKQzmZcvhz//hIMH9f0BAdC0Kbz8su445IizXbuKj4efftKZZr/9Bq6uuqh12bJQoEC6v7xhGNuVUinab5EgLMTTYNs2veS8ebNelnNzg8uXdaWlp0TtkWseuA9Z2MeTjR+9kOLnCQ0PZd6eeRyMPEjIqRBcnFz4/a3faVyq8e1rVs8mZNl4gv7Yh+kM0Lo1/PhjigLv5ct6efmPP3Ri1dWr+tdVr57ur9u0qePs7Tqcc+fg++/1LTJSH6n77TedcZ6BniQIy3K0EFmV1arfyceM0Zk6OXLoYg8JCfpd3Q4B+FT0KbK5ZSOPV54Mf21bZOSuOb6GhvMaYlVWANpXbE9ww2AKZr+j+e3MmZi6dMGUPTt07a2XPosVe+TzHj4Mv/+uf10bNuhfXYEC0KqVDrwvvigVqh5rzx4IDNT77s2a6Y3xoCC7Zz8/jgRhIbKqkBBo3lwHgHHjdGlJO7+TD1k9hA2nN3C63+kMf+3UZuQmWZIIOxuGyc/ElogttwKws+FMubzlKOjio5tVlCihM21feQVGjdJ1nXPe3yd52c4IvvnrMCf+c8cpvDBGeCEiTroC+nTSxx/reijVqmV4+9/MxWLRn1wuX9aFYypU0JVG2raFUqXsPboUkyAsRFZx/ryuyODpqd+M6tfXS3Evv3w7+9bOkpOy7OFJM3JjEmOYvmM6Y0PHcj7mPKf6niLIPwhPF0+93+vsStDqo9CqmK792LUrNGjAsogkRpmrcnbEhrsSjmJjYdiUy3w/R3H9SG2scW7gZCWbfxTvDnbhkx4+j5swC4Br1/TRom+/1YluVaroD5hOTrq4TCbjGP9nCiFSb98+PdOdN09n8rRvr+83DD2lchDnrp/jSNQR3qn8jl1eP6UZuSuPrOSbTd8QdjaM64nXeb7Y83zf9Hvye+enYPaCrO6wmpCZnxE0cx2m4zP1zHfAAKhX777kr9Nnzbz/eRQjrudm7xZP4uJ8cXJPwrPUBX0rHomTu5n9Pp4UK5byfemn1rx50KOH7ipRu7ZecXjttVQvOVuVlX+O/cOLJV60WxlVCcJCZGbBwfDRR3r226WL3gcrU8beo7rPwv0L6fxbZ5wMJxqWbGi3cTyqGIfFamFrxFZaLGxBvDkeZ8OZqa9M5d1q7+pqGGvXQs2amPxMmAq2hhdL6jO+AQG3nmPUqkNcj3Qj9nB+Yo8UIOFMblAGUTnjea8rLLy4GXe/KAznuxNiM7RSVGailM5nKFRILzGXKXN7v7d62jpw7Ti3g3ZL2nEg8gC/v/U7zco2s9Ggn4wEYSEyk8REWLBAvwE98ww0aqRnv927Q56MT3Z6lIORB3EynCjjW4ZK+SvRPKA5g2sPplL+SvYe2l22RWzjm03fkGBOwFTERJIl6dZjkTEXYP58ndy2c6eutNSli05wu8OBA7BkCYRNrEbiBb0P7Jr3GjlNR/EsfR73/Nf4Nrgp20fGEhF9/4mUDKkUlZkkJsIvv+gVnp07dQ/FiRP10a7581P9tFFxUVyIuUC5vOXwy+FHLs9czG85nyalmthw8E9GgrAQmUFUFPzwg34jOndOz35HjND7YQ5SdCDZtohtjNw4kqUHltK6fGt+fv1nAvIEMK/lPHsPDdBHjP49+S8ezh78ceQPQk6G4OPhQ8/qPXm+2PO4ObvpPV/lRFC/CbAz8nYxkzZtAD1B27kTFi/WwTf5/G72ouAVdACvMudxzRV76zUL3wyydqkUldmMGQOjR+sch3Ll9L/7du3S9JQnrpxg3OZxzNg5g0r5KxHaJZS82fI6RDtNCcJCOLqPPtLBNzZW95ibNUvPgB1MyMkQhq0bxpoTa/Dx8OGTup/Qq2Yvew/rLqHhoTSY24B4czwKRV6vvIxtNJauVbuS3T07XL2q93xP/kvQsHmYfArAn7qYiRUntmzR9R8WL9ZlIp2d9fndDz7QiejbLsYwZMmphwbZDK8UlVkcOqSXmg0DTp6EypX1Un/Dhmk6YrTr/C5GbBjBov8W4Ww407ZSW/o/199247YBCcJCOBqlICxMn3k0DH1otHVr/aZUybGWci1WC4Zh3OoSdDDyIKMbjqZbtW46qDmI6wnXmb5jOrsv7CbRkqg7G+FE75q96WfqB1u36hnYihWYTpzAVNcES3pj8fRm0yZY1E8H3ogIfcS6YUOdiNus2d27AIULPz7I2qpJRKZntepa5ePGwf/+p/fcn38exo9PUVGThz6tsmJVVlycXAgND2XV0VUMNA2kd83eFM7hgD93pVS63qpVq6aEECmQlKTUzz8rVaOGUqDU6tX2HtFDxSXFqSnbpqhS35ZSfxz6Qyml1LX4ayo+Kd7OI7vb+evn1SerP1E+I30UX6Ba/txSeQ73VM5fOivP4Z5q07wRStWpo3/eOXMqNXiwMl+IVP/+q1TPnkoVKKAfcndXqnlzpebNUyo62t7fVSYXH6/U5MlKlS2rf7iFCin19ddKRUam6WnjkuLUtO3TVMCkAPX9tu9v3Xc1/qotRv1EgDCVwhgpM2Eh7C0uDqZM0c0UTp3SWaDffQc1a9p7ZPe5Gn+VKWFTGLd5HBduXKB6oeq6OxA4xMw3NDyUkJMhBPkHsebEGoatG0aiJZGW5VoyqNYgahapefsawx9T3Tbg74917Hg2BXTmpz+ys7giXLyoE86bNoXXX9dHraViVRolJuplBKVg6FBdSnL+fP0DTkPx66i4KL7f9j0Tt07kwo0LVC5QGb8cuouFh4sHHi4eNvoG0ofUjhbCXhISwN1d/1ms2O1mCq+8kqbluPSilKLyD5XZc2EPjUo24sPaH1Lfv75duiE9yMjVv/PxhtdRyoKT4UrL0t3xzRHPANMASvuW1glt332nKyx9/z1WK+z7bi0zDtZm4RIXzp8HLy/943/jDd3T4gkaHImH2bpVLzHv2AH79+t/2xER+tiRDf7tNJjbgDUn1tCkVBMGmgbyQvEX7P5vUmpHC+HIwsL0/uO2bTqt1t1d173Nl8/eI7vP0aijTAmbwv+98H+4u7gzssFI8nvnp2rBqvYe2i1KKYau+pGRmwejjCQwwKqSWHvwKlNbDKf0uSsw8B346SdUUhKXg1rxf70tLFzszNmz9fDw0DPe1q31nxJ4bcBshmXL9H7vpk26bnnXrnrVx9sbCqd+bzbsbBjjN49nXONx5M2WlxENRuDp4knF/BVt+A1kHAnCQmQEi0X3oktuppA9uz5rGhen/9vBAvCOczsI3hjMov8W4eLkwmtlX6Nusbq8VPolew/tFqUUP+75kVGbRrHv4j6cVE7AGVAYuOBsLs+B4eNovmQMFg8vtpR/l4/O92X9v6Vw36Rnuq1b65mvLDXb2KpVejmhRAm9zfLOO2n6IVuVlRVHVjBq0yjWnlpLDvccdHy2Iw1LNqRG4Ro2HHjGkyAsREb43/9uN1MYO1YXfMiRw96juk90fDRvLnqTv4/9TQ73HAyqNYg+Nfvc3SXIjkLDQ1lzYg0vFH8Bk5+JqdunApAnsT9eluexWg/gH/UHbtZKHHQOZO5pZzzy5GZYZHeu781No0YwZ6T+VTjgjz/zOnpU13IuXBg+/BCaNNFNkBs3TvPWSmxSLDWm1WD/pf345fBjTKMxdK3alRzuWeMXKEFYiPRw7pxuppAjh35TathQL881beowzRSSWawW/rv0HxXzVySne05cnFwY2WAk3QO7k9Pj/i5A9vLn4T9p8UsLzFYzHus9WNNhDUvfXEoerzy8PHQZ9TcsoWPYn+SPjWKmR0m6xAdx1lD8Xq8uI9+Gli0drqhY5pZcynPcON2D0cUFevbUjzk762y2VLoSd4U1J9bQ6plWeLl68VKplxhSZwity7fG1dnVRt+AY3CsdwMhMrs9e/RM96ef9L5Y5876ficnXWjegSSYE24t5569fpZTfU+R2zM3y9ssT5fXW7YzIlVFKg5FHmJM6Bhm7pyJRekiGEmWJEJOhmDyMxH78XCWjR6Be1Isq2jEGAaw1ieQfBUPMGKgD50bOcYsPssZNEhvr/j66q5dPXpAwbT9rE9Gn2T85vFM3zGdOHMc4f3CKZS9EKMajbLRoB2PBGEhbGX4cPjsM51i+957ush8yZL2HtV9ridc54ftPzA2dCznYs5RtWBVulccxSvjt3PuamK6VHG6t7tQRHQcQ5bsBXjk65yPOU/5yeVxdXbl1bKvsuLoCpIsSbgZLlgP1eLVUVB8uQ+Vra2ZX6gXh8v7YC1xkkD/XTe/B9sG4NR+kMgSIiP1UbrWrXUm/1tvQdmyuqSkZ9pqX5+KPsWH//uQX//7FSfDibcrvM3AWgMplL2QjQbvuCQIC5Fa8fH6nGOdOvrNqEkTvST33nuQK5e9R3cfpRSGYXD8ynEG/TOIBsUbMLfFXK5HB/Dx0n3EJSUCKQ+QT2LUqkN3lXIEiEuyMGrVobteY+Ppjfyw/QdcnFyY+dpMCngXYE7zOTQs2RBfl9zMHx7Mnm1TaLX/DFPOnGRH4XqU7f8BldtAp8rJJ15K2GTM90rtB4lMb/9+fcRo3jz9bz5nTh2EAwP1LZWUUlyOu0werzy4u7iz5sQa+j/Xnz7P9aFIjiI2/AYcm5wTFuJJXbwI33+vz5xeugRffOHQzcSPRR1j9KbRJFmTmP7qdAAOXz5MGV/d8rD2yDVEPKCVXmEfTzZ+ZJset8U/Ws6D3mkM4MTIpiRaEvlq7Vd8vf5rFAoDg/+1/x/1i7/A7u1mTg6aRNX1EyhqOckxp1Jsfq4fRT7pSJ3G2TLsSHVG/JwcilJ6I33ZMvDwgA4d9OpOuXJpetoEcwIL9i1g9KbR5PbMzbp31gHophnOqS/a4UjknLAQ6aVfPx2AExJ0ktWAARAUZO9RPdDOczsJ3hjMr//9iouTC12qdLk1G04OwPDwXra27HFbyMfzgQGskI8nG05v4M1Fb3L2+tlb9zsZToyZt57eC19g/35ndjGbGN8ihHYZR5WhzWibLeOLmWTEz8nu4uJ0VvPrr+tlhQoVdNvMbt3SnNUWHR/ND2E/MGHLBM7FnKNS/kp0q9bt1r/JrBKAn5STvQcghENTCkJD9Z+gsz47dtQNZP/8E+rXt0nVH1v7IewHqk6tyoqjKxhUaxAn+5xkctPJD6wk9LBetrbscTuocVk8XW8HTjNROLmGM6hx2Vv9hkfUG4ur4YlhdcY5EXpPG0Mhr2gmTzbwO76WZyLXYwpujocdAjBkzM/Jbs6e1clVfn76fO/Onfr+YcPg449tklY+d/dcPlr9ERXyVWBVu1Xsem8X7Sq1s3t1K3uTmbAQD5KQAD//rDOd9+zRBTbq1NF9Th2QxWph2cFlFMpeCJOfiVfKvMKI+BG8H/j+Y48ZZUSP2+Q90w+XT+Fk4s8kOp+gXO5qNKvUnX//hbwrlrN10QrG5C1HjP8OakV48uyLXfh7lBnyANj/qFSW7AV84QIMHAi//KKz+Zs310vOlSun+al3ntvJ6NDRNCrRiI6VO9K5SmeeL/Y8lQuk/bmzEgnCQtzpxg1d4WfSJH3Wt3x5mDEjTQko6SneHM+Pu/UxoyNRR+jwbAdMfiYK5yjMR3U+StFzZESP27CzYYzbPYDDlnXgDM6GC2XPD8TfH86cgZreB9ic8AoJiUVwa/YNRrd3wcfHZq9vC1mmF7DFon/oxYrpKlYbNujjRb176wpXaaCUYtWxVYzeNJrVJ1bj7eZNYEH9/463m7cE4AeQxCwhAGJj9dGiuDj95lSlCvTvD40aOeRyM8C07dMYGjKU8zHnqVawGh/W/pCW5Vri7OQYzR/u7Os6aeskBv0zmHhzPKDA4swLIUH0Ol2BhODxvPoqeG74R++vu2atYgwO4/p1mDVLf8h0dtZ1y52c9AzYRgVk2i5py097f6JQ9kL0qdmHbtW64ePhWB+mMoIkZgmREkrBunW64MDBg3qf19NT/3fu3PYe3QOdu36OPF55cHV2JTo+mor5KjKvxTyH6ByT3CKwll8tjl85zqhNo+hdoy8lrnRj3eyumMOewWjdFMM5AXerheEnVmMKzA5vWHUwaNjQruPPssLDdeCdNg2uXYNatXSCYfIELA0B+Gr8VabtmEaXKl3I5ZmL9pXa07BEQ9pUbPPUJlo9KQnC4umTlAQLF+r93h07dMWfHj10v1NPT4cMwEcuH2HUplHM2T2H6c2m0/7Z9gysNZBBtQfZe2iADsAN5jYg3hyPunkYKZ96lk/7FOZyKOTK5cHsyicoMSeef0u5UL98U0z/jISAADuPPAtLnuFu3qzP+b7xhg6+NdLe8CD8ajgTtkxg6vapXE+8Tv5s+Wn/bHualGpig4E/XSQIi6fPX3/pKj8BAfDDD9C+fZor/qSXbRHbCN4YzJIDS3BzdqNz5c7ULlobwO4z3zuFnAwhznzzqI4Cdr3Dld+nMLzir1QbFkrtgSY8LjaEmZ9j6tkT8ua163izLLMZFi/W9ZybNtUV3Fq0gBMndOZzCj2sMpjZaqbzb51ZsG8BSineKP8Gg2oNcqjWlpmNBGGR9R07ppfjihSBwYN177qVK/Xyp5PjntJTStH5986cuXaGIXWG0Ltmb/J757f3sG45cOkA4zaPp1XOrwn9JQhyu4OTGcPqxgRXL94vUBKXPWegRlfwMEHRorqwibC96GiYPh0mToTTp3W51GLF9GMuLk8cgO/MAj8THUv/JX8AzWhepTBx5jg+qP4BfZ7rg7+Pv+2/l6eMJGaJrEkp2LhRLzkvW6bfiHr10vu/DspsNfPr/l/5YfsP/P727+Rwz8F/l/7DL4cf2d0f3os1o+oZJ+/55vbMzbL9f7Hy5O8YZk/Uz4vxOvsSQe1DeVak8P0AACAASURBVCbbF7T4Yz21jsTpM9QDBujGvQ78YSdLeOMNWLQI6tXTS86vvJLqFoLJlcEUZm44r+eay1KSjJNUdZ/H9iFv3yquIR5OErOE+PBDGDVK7+8OGQIffJDmDi/pJTYpllk7ZzEmdAwnok8QkCeAU9GnqJi/Is/kfeaRX5tR9YyT93xvLTnHZ4fNn1PZ3JPPWp8i/sVTTNoSh8/fRbierRb//jSI+m83ttnrizsopc+tjx+vz62XKAFDh+qiGlWqpPnpz0RHcd1lJdecf8fidAlXqx++SR9wOV5v2UgAti0JwiJrSF6Oa94cSpXSZfeKF9fVrby87D26h7p04xLlJ5fnUuwlnivyHGMbj+XVsq/iZKRs5pjSxgiplWBOYOmODUxfsZU4a6KusWd14jlrP35uXplii1rC1xsYtH8QEQH1mFajJdMAz/8UI3ZGZL4ztI4sMVEnFI4bdzuh8OBBHYQrVkzz01uVFSfDibw5LZxOmIO7NYDcCT3wtFbDwInCWaEymANKURA2DCM/sEgpVdcwDFdgCZAbmKGUmpmeAxTikU6c0Pu9M2ZATIxedu7bV2eA2iALND2cvnqajac38nbFt8mbLS/dqnWjUclG1C1a94lnGelVzzgyJpoBC37gl5MTSHC5AAt/xekNNyARdycnxq6bQbHNEVCsGBOa9uCv4tXv+npbfhAQ6AAcEKD/vZcrpxMK27WzyQfM/Rf3Mzp0NBdvXGR5m+V82qQOA5ZMx5J0u1Rlpq8M5sAeG4QNw8gFzAGy3byrF7BdKfWFYRh/GYbxq1LqenoOUoj7KKWzmhcs0PuNb7+t98JssByXXvZd3Mc3G79hwb4FuDm78XLpl8npkZPhLwxP9XM+qjHCkwoND2X+tuWE7D7Bf+Y/UG7XcbvQkLcLzCVoWEkmHQ4m/EYY0//4m7IqL/wyFlq2ZPynqx7YISlLNTbIQMl7/B7HjtD87E78hn+mP8z06qUDcKNGad5jV0qx9tRaRm0axV9H/sLL1YvOlTtjtppvfnB6MfNXBsskUjITtgBvAr/d/HsQkFwPbx0QCPxr85EJcS+zGVavhsaNdRWrggX13m/PnlDYcd8gDkUeYsDfA1h+ZDlerl70rN6T/qb+j63pnBK2qGdsNsPIBev4/EgTrEYiGBZ84xrySYVger3mSsTQYWQb/zcT3p1KTvfWDKvfhE+z+zCidCWau7jY9IPA027ZjjMsGzuPrzYvocGxbSQ4u9KkZG3gBZr362ez15m9azadf+9Mvmz5GFZ/GO8Hvo+vl++tx5tXKSxBN4M8Nggrpa7BXZvx2YCIm/8dBdx3ZsIwjG5AN4CiRYvaYpziaXb1qt7v/fZbffxi2zZdy3nUKHuP7KGsysqVuCv4evni7uLO9nPb+TLoS3pW73nXm11apaWe8alTis9nbOTnM9+QkH0/+CSCkwVnw5kBFYvQb/ZH8Pff5Hd1Z2GFF3GzJHEDiPbMAWbrreXmLNnYwB4OHKB801dofv44kV45GV/7beZVeZlIr1xpXtqPSYxh5s6ZFMtZjNcCXqNluZYkWZNoX6k9nq7yYcmeUpOYFQN4AlcB75t/v4tSaiowFfQRpbQMUDzFoqJ0K7Xp0/V+b1CQbqxQ1XELAySYE5i/dz6jNo2iuE9x/mr7F/4+/oT3C8fFKX3yIJ9k1mKxwF8r/p+984yK6urC8HMZehNsIIjYe0NRihoxauy9RWM39t5N/0zsKPZYYtfYYtckNhRNADvW2BsqihUR6TP3+3GCJsZCGRhGz7NWlgoz924mrNlzdnlfHWPXbeW45RTIfwgzl1zUdW7On09XC2N1E1P8vl0Kyc4wbhw+9wvzxMr+P9dKKTe/N8YGhuDBA2Gm4OEB+fMTpbFkZINBbCvtR4LpS9nH9Jb2I2MimX1kNj8e/ZEn8U/oXrE7zUo2I4dlDnpV7qWvn0KSAdLzrnAcqA5sACoAh/QakeTDRlVF8s2VCywsYPVqMfE8dGi2Tr7RCdEsOLaAGYdnEPEsgorOFelcofOL72dWAk4tO06GMnNbEGe2+RHJGWjSGwe1EMO85jLcoynWi1cQGtOSoIZl8Cvoh0+Fh6L3aGGB9aR9PHlHuVmWL9PIX3+JFaOVK8U0/+nTYGfHkAGz9Vban/jHRMYeGEuiNpHmJZszwncEvm6++oheokfS886wHPhNUZQaQGngsH5DknyQpMjtBQQIkflz58DGRkyDZuMVoxTmHpnLl/u+pHah2ixttpS6hesafJ9SVeHAAfhu+W4OujUGRYtJQwu+KbyNkh7raGtdAdOZs6HVSIiNxadJE3zGjxZDP/8QWJLlZj0SGgrffy8U2ywtoXNnMc3/9+9KRl5rVVUJvhVMeafy2FvYU9ChIN0qdmOoz1CK5yqeaT+SJGOkOgmrqur39583FUWpizgNf6uqqvatT5RI3sar/d6iRcWpV6sVySCbJuBLjy4xNWQq9YrUo1XpVvTx7EPdInXxdDG87/DTp+KANXPpba7kngFVZ4MmCQBFScSqyFE6HM0DvUqJla7PPhO2jW/YNZXl5gwSHw86nfhdvn4dwsJEm6V37/9oaKfntdbqtGy5sIWpoVM5dPsQ0+tNZ4j3ENqXa0/7cu0z9UeTZBwpWykxDKoqPv2vWweffir6vcOGCdH5bCxxePj2YaaETGHz+c2Ya8z5odYPGXYy0pfs5MmTMG8e/PwzPK/8A4rf9ygmKn7uHxN8+yDJ2iTMTcwI7LofnyQnsVudjZXEjJ7798X/kB9/FL/bo0eLio9WK1otGURVVeYfm0/AoQCuPL5CYcfCDPcZTteKXbE2y54fXj8UpGylJHuiqhASIhR/qlQRb0qtWgn1n2y835tCj609WHJyCQ6WDnozVMio7GRioqjiT1wVwhmTZZhdaU3Htp9QvHkp7pr1Y2i5nhTcGEjotjMEWdzDr/wn+Lj5iCePH5+h2CVv4OxZ0e9dtQoSEsQHy2rC+QpT0wz594KQObU2s0ZRFDac30BOq5z80uYXWpRsgcYkfXrREsMhk7Ak80lKemmvduQIODqCt7f4nqlptk3ASdok1p9bT7OSzbA1t6V+0fqUyVuGnpV6vtVQIS2kV3YyIkKIJs1fqOV+0SlQ+2tQdKhVl9Gz6wF83FrDlGvQvjo8fYpPtWr4DPsRmjbVS9xpJUmbxC9//cLWi1uZXm86LnYuab5GVhlVZJjhw4W2c7duot9bQj+98yuPrxAQGsDqM6s51+8crvaubG63GTtzO4PPH0jSj0zCksynZ09Yvlz0e+fOFXrONjbvfp6BiEmMYdGJRQSEBnAr+hZLmi6hm0c32pRpo/d7pUV2MsUYas4c2LBRRVthEdZd/cHy8svH6LQEXdsnTrtarZhwHj4cvLz0Hnta6LOjD0tOCoVbb1dvhvqkTXgiq4wq0kxcnKj/z54N27cLu8a5c8UHzVz62Qc/fPsw/iH+bDq/CTONGZ3Kd0L9W6PM3uK/q2MS4yL7Nt8kxsvVqzBokBhCAaFotW0bXLwI/fpl2wScrEvmm33fUGB6AYbuGkohx0LsaL+DLhW7pPuaW8LuUG3SPgqN+ZVqk/axJezOv77/ptWTf349NlbMrnl4QI1a8ezaBUMGK1Tr9QulCtnzQ82xWCnmaHRgnqTD78bfTxwzRgj+GyAB34m+w+g9o7ny+AoAA6oOYHv77RRxLELQzaA0X+9tFQODEBkJ330nkm7PnmKO4f598b2iRfWWgO8+u0u1JdUIvB7ImOpjuDH4BouaLiK/fX69XF9ieORJWKIfVBX+/FOsGG3dKsrMXl7CyahKlXc/34A8iXuCo5UjpiamHAw/iF9BP0ZVG4V3fu8MXTc1p7e3raSEh4uZnp9+gsd5tmBddyoWLU5zqPdZSjgX4GnMGuxXb0TpN4PazxIJKm+PX41O+DTrLy5kgBLl6cjTTAudxuozq9GpOornKk7RnEXxyOeBRz4Ptl7YyobzG9DqtGnqX2aWUUW6iI6GIkXEp6PGjcXQVc2aenm945PjWXV6FSfvnWROwznks8vHtvbbqFGght5aIJLshUzCkoyj1cJHH4mhqxT/3v79wSXtfb+s5FjEMfxD/NlxaQdXB13F2daZPZ32YK4xf/eTU0Fq+r2vrqTky2FFozzlWDU+D602g5r3DE5dR4P978QCGkXDufuhlHAuQA5rR5g2DWxs8PlqFT5t24KZmV5iTyuqqtJ8XXO2XdyGjZkN/av0Z7DXYAo5FvrX40b4jmCw9+BUWzWmYFB9alWFXbvE0vXEiWBvL3oCvr5QXD/7t0/injDv2DxmHZ5F5PNIPJw9eJ74HBtzGxoWa6iXe0iyJzIJS9LHkyfw66/CTk2jgdq1hatR587ZdrcXRLLYfXU3U0KmsO/6Puwt7BlQZQBmJiJ56SsBQ+pPb809XKlX0pU1a8S69NenxGeZ/iMfMt+mMo8ARaeIPqBOx8UvesH6xqKsf/Ag5M1rkFNvkjaJPdf20LBYQxRFoUyeMni5etHHsw85rXK+9jklcqdvSMkggiFxcWLCefp0OH9efKgcOVL8z+naVW+32XN1Dy3WteB50nM+KfIJI31HUrtQbTls9YEgk7AkbVy+LPx7ly4V5biqVcVp4PvvDR1ZqrgRdYMGPzcgn10+/Ov606tyr7cOt2RkIjc1p7eICFFynr0llOic+3AqqlCn40O29gvA2jo3H59fh+X5S7Q8+SWJqJjrVPwK+YnX3sYGnDK2IpUeohOiWXh8ITMPz+R29G2O9TxGZZfKTKg9IVXP//XSr1yPus6AqgNSfc8sFwwJDhZyqQ8fimb8ypXQti2Y6+dDWtjdMOKS4/B186WyS2XalmnLYK/BVHCuoJfrS4wHKdYhSR3h4cLPdPt2UfLs0EGsX1TI3m8azxOfsyRsCRcfXWROwzkABF4LpIZ7jXeeel/t6YI4fU1sWS5Vb/5ve35+rSszZgitkmTXIJQun6CaCFUrVztXzvU7J6wOw8KgUiV2FbNhdtUiPC7UjlEtu6Qp+ehrtScqPorxB8ez8MRCohOiqVWwFiN8R1C/aP00lZd7be/F+nPreTTqUfbaaz1zRvR7q1WDqCgxcNW/v976vaqqsuvqLqaGTCXweiA13WsS1DUo43FLsh1SrEOiHxIT4dYtMYTi4CBECL7+Wkw4OzsbOrq38uD5A+YcmcOco3N4HPeYGgVqkJCcgIWpBbUL107VNdK7w5vCf/q9dtZUNy/P1IG5CA4GOztoOTiUHY6NiU0WCdgEE/rGlyPHxAAYO5Yt5GVfq6/Y5V6JBDMLSCJNqzn6WO15lvAMOws7TE1MWX5qOQ2LNWSEzwgqu1RO1fNfxa+gHz+d+IlTkaeolM/Aphwp/d6AANizB3x8xGyDgwP88ovebrPj0g6+DPySM/fP4GLnwuQ6k6WLkQSQSVjyOh49EkoQc+eK/tfp02IY5fLlbC0pmcKOSzto80sb4pPjaV6yOSN9R6bLPUYfE7nNPVzxK+TK4sVildT/JhQoF86wqff5rqcnOvNSNFtbmcO3DglZyWQdHy/aCeVNQFXx33WRO0V9/nXNtHwQSO8HCVVV2XttL1NDpxL+NJxz/c5ha27L9cHXsTHP2IqZX0E/APZf32/YJLx1K3z5pXA0cnERQ1e99JcYnyU8w9TEFCszK+7F3EOn6ljabCkdynXQ6+yBxLiRSVjyksuXxbTtihViKKVuXbF+kUI2TsDHI46TrEvGK78XXq5edCrfiWE+wyiZu2S6r5nRidzr10X7fOFvocQ5BVHSKz81h+0m+OlaDthVwN7+GODAAW1nQhf9SVBBBb9idfHZOvGFbWNGPwik9flJ2iTWnVvH1JCpnIo8hbOtM4OqDiJJm4SFqUWGEzCAi50LxXMVJ+hmEMN9h2f4emkiMlL00m1txXChubn4fW/XTm/93jvRd5h1eBYLji/g+1rfM8hrEF0rdqWHRw85bCX5DzIJf+ioqhCVNzODQ4dg2TIx8TxkCJQta+jo3oqqquy5tocpwVMIvB5I3cJ12d1pN3ls8rCwycJ3Pv9dvdL0TuSGhorq5qZNoLiFonapBSYJXABuRFsysEp/hkSXEY4LFStCjRr4tB2Gz6BB4Ob2r2tl9INAWp+/9eJWOm3uROk8pVnSdAkdynXAwjTjZgOv4ufux5GII6iqmjWJ6dw5MeW8apU48Q4dKqb5u3TR22T52ftnmRoyldVnVqNVtbQu3ZoaBWoAhveTlmRf5GDWh0pCAqxeLd6YOnUSqxeJiWIgJW9eQ0f3TnZc2sHX+77mVOQpXOxcGOo99J2Tzv8ktUNXqR1q0mph82ZRSDh0CHI4aunZOxnT6gFMPvYVKioKCl9b1+f7hZfhyhUx+LPw7R8WMnM4rLmHK3ei7zDz8Ezy2+dnkNcgknXJBF4LpG6Rumne5U0LKf35TGf3bvGJaNcusLJ6qedcrJjeb/XR0o84fvc43St2Z6jPUAo7Ftb7PSTGgRzMkryZ+/dh/nzR771/X3jIFikivmdunq0TcExiDGYmZliYWhD+NJwkXVK6e2yp7ZU293B9a7KLiYElS4RpzvXrUKh4HO38l3PUbCouVfvjnd+PmSctSUyOxzxJpcGs38GlqhiLbtnynXFmdDXnTc8v4vKYLlu+fKFs1b+KUNkyNTGlXtF6qbp2RsjUBKzVit11gEmT4MIFmDBB9Hv1JCeZrEtm8/nNzDoyi1/a/IKzrTMLmywkj3Ueclnr5x6SDwN5Ev7QaNQIfvsNGjYUJbnatQ0i9JAW7j+/z+zDs5l7dC4Tak+gj2cfknXJmCgm6T6tFRrzK6/7zVeA65MavfP5d++KQat58yDKNpT8tX6jlGckJxO28CD2AVVdqzK2cHfq1/yc0IgjBC36Gr8ryfj0HS9WYAz4mv8v6H+MPTAWazNrenj0YKj30P8oW2UFI3eP5FHcI5Y0W6KfC6b49y5aJMoRrq5w+7b4YKmnfu/zxOcsO7mMgEMBXHtyjWI5i7Gq5SqqulbVy/Ul7wfyJCwR6HSiDDdrlih7urmJfti0aVAy/QNLWcWVx1eYFjKNpSeXkqhNpHnJ5lRxETrUGe2xpbfXev48TJ0qWotJSVCzYyihxWpzWxfH7Sfgm9+XDYVGU2PRbpRdfWCzEz7Nm+Pzv70GS7xJ2iQ2/LUBXzdf3B3cqVO4DmYmZvSt0veNylZZwbPEZ2w8v5GFTRZm7P/nX3+JUsSKFaLN0rAhPH8uvpdff0YHzxKeUWRWER7EPsA7vzdT606laYmm2WvXWWJ0yCT8PhIbKxR+ZswQpTgXF9GDdHOD8uUNHV2q6by5M8fvHqdLhS4M9xmebsnD15GWoStVFeqQU6fCjh2itdi8zykSK0+jpKsrfwQnAqBBofHOa3y0YYTYox43DmqIwRxDJOBnCc9YHLaY6YemE/40nLF+Y/m25rdUL1Cd6gWqZ3k8r+JX0I8Fxxdw8t5JPF1SdWj4L/fvC8EYU1MhJTlkiF4/YF5+dJm91/bSt0pf7CzsGOE7gmpu1ahWoJre7iH5sJFJ+H0jNhYKFxarGJUqiSNbmzZ6K8dlFqqq8vuV35l1eBarWq4it3Vu5jeeT16bvDjb6l8YJDW91pRhqylT4OhRyJ1Hpcv/9hNeYArrw3dhe8eWaoW+xlxjTqI2UdgI3rcWTeIOHcAiCwaP3sB3+79j1pFZRMVHUaNADeY0mEOj4u8us2clNd1rAhB0Iyj1STghAdasEZPlM2aIUvPq1eDnB3ny6C220Fuh+If4s+XCFixNLWldujV5bPIwqtoovd1DIgHZE34/OHFCqP2MHi3+PX06VK4sTmHZvN+bqE1k7dm1+If4c/b+WfLb52d96/XClN5AxMfD8uXi5HvlCrh6h1KpVSCXLVdz4dF5nGycGFysI30Cn2K+aTc1eg3meuIpyumKMqx1N5pXMozX642oGxR0KAhA963deZb4jBE+I/DKn/V+wqml5JySFM1ZlB0ddrz9gQ8fioHCOXPEB8xy5UTfV89mIecfnKfXjl78Gf4njpaO9KvSj4FVB+Jkm/Ua3RLjRfaEPwS0WqHjPH26qJXa2UH37uI0MHSooaNLFc8SnlHmxzLcir5FubzlWNF8BZ+W/RQzjWHs+KKixFzPzJnifd6jaiwd581l46Pv+C02ESVOYbR7R/634xmWowLQaUzZUtoPHjqTw7oE4cAXm8+ComSescArqKrKH+F/vLBkPN7rOJXyVWJR00WZumKkL7pU6MLzpOdvf9Dvv0OrVkJApn598ftdt67ePmDGJ8dzL+YeBR0Kkts6Nw9jHzKj3gx6VOqBrbmtXu4hkbwJeRI2Rk6cEI4uV69CgQIwaBD06CH0brM5Ec8i2Hd9Hx3LdwRE2dQ7vzf1i9Y3mJrQ7duisrlggVg5qtXoEe6t57LjwWwexj7ERDFBp+rQoOGHvVq+OOcIffvSVK3Iad1/T2KuDlYEj/k4U2PW6rRsOr8J/xB/jkYcJbd1bgZUGcCAqgOMf0VGVSEoSCi01awpZFS/+koYiJQpo7fb/NPDt5BjIUK6h6AoStYJiEjeW+RJ+H3k5k148AA8PUXP191dTDq3aCGGUrI55x+cZ2rIVFaeXgnAJ0U+Ia9NXsbWGptp93yX0MbFi6Lfu3KlGCRv1uEB5nV+YNvtxey/GUvjwvVpeNeO4XFbSFQUzDXm+HUdDZ+OABsbzoz59bX3TYu2dHp5nvScz7d/Tl6bvMxrNI8uFbpgZZYFBveZQLIumUexj3Ayd4T164W4RlgY1KsnknCuXKIUrSduRt1k+qHpLDqx6F8evinIBCzJSrL/u/eHTmioKDlv2iT6vIcPixNvYKChI0sVVx5fYdiuYWy/tB0rUyt6VurJMJ9h5LXJXFGQt7kHuelcmThRvKRmhUKpOHg3o9t8gl/Z4hSdvZK2BRsx4rQdZQZuhKdPqdiyKkFDmuFXsNa/etUZlZRMCymuUH+E/0Fg50DsLewJ6R5CydwljX5FxnuRN05RSfw686EwWC5VCn76CT77TK/3STnh/nr5V+YenUv7su0Z7jNcevhKDIpMwtmV3bvh229F0s2RQxgpDEi9Cboh0eq0PIx9iJOtEzZmNhy/e5zvan5H/yr9yWOjvwnWt/GqIpaqwuMrDnRta8XTK2CfQ6XOV7MJNBvKUVVHpz2TCcwXyK2c47DtMkQcjVu3hmHD8PHy4nVjYunVlk4Llx9dZlroNJafWk58cjxNSzTlacJTHCwdKJNXf6VZg3D1Kjg5UcWlCj/fW0ZyaV9MFy0SJ2A9mYWk6Iv7h/jTqlQr+nj2oWvFrjQp3gS3HG7vvoBEksnIJJydSHF1sbERPr6PHolp0C5dhOtLNic+OZ4Vp1YwNWQq+ezycaDrAfLZ5SN8SHiWn9ZSSsKqCnGXnXh6qAiJdx1RbGNoP34t5xynsOd+GCmyWYnJiQTdCMLHt6n4sDN4MBQs+NZ7ZFRS8l3sv76f2itqY6Yxo3P5zgz3HZ4hV6hsgapCcLAoOW/ZArNm4efnx/zj8wlbNokqrlX0cpskbRLrz63HP8T/hRtU+7LtAbA2s8Y6h36nqiWS9CKTcHbg8mUxkrt0qdC4HTxYJN5u3bK1fWAKj+MeM+/oPGYdmcX95/fxdPGkf5X+L8p/hiiX5rOz5tIhR6KuP0TruAdNjmfkLJcbx2o/syZpMsWTizHGtgEzo3aSqKiYK3/73LqVEeX/VPIubem0oFN17Li0g/jkeNqWaUu1AtX4vtb3fF7p80zZlc5SVPVlv/fIEXB0hC++gJYt8bMXvx/7b+zXWxJut6Edmy9splTuUixuupjPyn2WNYYREklaUVU1U/+rXLmyKnkD+/apapMmqqooqmpurqpduqjq6dOGjirNzAidofI/1AarGqj7r+9XdTqdwWKJj1fVhQtV1Sl/kkr+EJWvLVW+U1S+M1Xdvw5Q1x29om6Z0l3V5s2jqqCG1CqmTpjVRg25esBwMSfFq4uOL1JLzimp8j9U38W+BotF7yQkvPx79eqqWrSoqs6dq6oxMf96WKk5pdQGqxqk+zZ3n91Vvwr8Sn3w/IGqqqq6//p+dfvF7apWp033NSWS9AIcU1OZI+VJOKvR6V6ebseOFT6nX38N/foJqUMjIOxuGP4h/tQrUo8uFbvQo1IPahWqRXknw0lixsaKWR5/f7hzB8rVuIWm7iAidPHiAWoyNcreo61nEdiYFzyrwPDh+NSqhY8Bp2HXn1vP4J2DuRdzDw9nD1a3XE2bMm0MFo/eCA8X1Z2VK+HMGXBygg0bxB77a6o7k+tMxsEy7St2Fx5eeDF1n6RNolzecrQr205UNSQSI0Am4awixUJw8WKh9JMvn5BlyptXiBFnc9R/DLjsvbYXO3O7F84xtua2mZKAU+PlGx0tXBmnTxcbXB99BHW+H8fK29+hqAoaTEDVYZ4M/ZS/J7LHjzdomT/8aTiWppbktclLDosclHcqz8oWK6ldqLbxr8ccPSoMQjZsEP9u105ITYJIxG+gSYkmabpNkjaJthvavpCV7OHRg2E+wyias2h6I5dIDIJMwpnN2bNCCWLVKvFm1KABPHsmkrC7u6GjSzWdt3Rm1elV5LPNx6Tak+jt2TtdJ5fU8rYVo+Yerjx+LA5aAb+EEpNrP6Ua2bC4U3uafJyX3y9UJO/uBgxed5Pw8LMElbHGz/tTfPw6iYsbKAGfuncK/xB/1p5dy1Dvofh/4k+9ovWyxL83S7h9G7y8wN5eTPMPHChMQ1LJ2N1r+eXYPZ5HF3u9lrdOy9GIo3jn98ZMY4ajpSPffvQtA6oOyLKpe4lE30jFrMzk3j1hpWZuLgatBg82CgtBEJKSS8KW0LlCZxytHPn98u/cjbmbZQMu1Sbte+0Obl5Te+roajBnDjxz/AOlax1UE+Fi1MujJwuaLhQeg0WKCF3hYcOgUyeDMHo+8wAAIABJREFUVhv2Xd/H5ODJ7L66G1tzW3pW6skQ7yEUyFHAYDHphefPxTDhpUvCLhNg61b4+GMho5oGtoTdod1WX0x0eXBKFAIuVmYaJrYsR72yOVlxagXTQqdx9clVLg24RJGcRfT900gkeiNTFbMURXEEfgbyAsdVVe2d1mu8tzx/LjxN//pLOL47O8O6dVCrFuQ0nG9rWrj77C6zj8xm3rF5RMVH4WDpQJeKXWhQrEGWxvGq6pQ2xoLoI4UJP1mAY8ng0ftHzrt8SZxOJGBFhTyLN0C92cK96OBBIelpoFOvTtW90G5edGIRpyNPM7H2RHpX7o2jlaNBYtIbERFidW7+fLFW5+MjqjwWFtCsWbou6b/rIubacjzX7EMlGQVTnidFM+jXb4jfu4MHsQ/wdPFkbau1L0wqJJL3gfSUozsBP6uq+rOiKKsVRfFUVfUDPer+za1bojG5cKF4U/LyElY8lpZCeN4ISNYl03dHX1acXkGSNomWpVoy0nekwRx4UtSokqMtiT5SmJhTBVBNYslT9gEHV+ej64FtmIVbk6SJRkXFXAf2UYX4LfgiDT8u/84d38wiJjGGxScWM+PwDLZ+upXyTuWZUX8GOSxyvB8rMps3iz5vcrKQTB0+HHx9M3zZiKg4LE3KEWP6Kwkml7DUlUZHAneSV9GwcF1G+o6kpntN4++ZSySvkJ4k/AgoqyiKA+AG3NJvSEbG5s3Cr1dVoWVLYSru65vtLQRBDFtdfnyZ4rmKY2piysO4h3obcEnNUNXb6FahNCO/SSDqwS0oORdN17/QuQTxXZ1dlCyZj1Kz6nFo3i72FrJgrndRLudpzo8+PrgeeUjDzPVOeC2RMZHMPjKbH4/+yJP4J1QvUJ1ErTilZ7ZEZ6aiqrBzpxCLqVEDqlUTk/wDB4qSv55wcbAiPKosAJEWo3CP24EpOfG0WM2vHVrr7T4SSXYjPUn4T6ARMAg4DzzWa0TZneRkITqcKxfUri3emIYMESpLBjp9pRWtTsvWi1uFA8+do1wddBV3B3c2td2kl5PGu4aq3kZ4uPClWLzYGW25JdC1Jyg6tAp8pC1J40M7oUYNguyKMaLhEH4v7stzi5fqR1lhnvAqCckJlPmxDI/jHtO8ZHNG+o40qB+yXoiPF8OE06eL9krz5uJ3PW9eMWioZ4QEaCL2SS0BUNFibWbOV/WN/HWUSN5BepLwd0AfVVWjFUUZBnQDFv7zAYqi9AJ6ARQoYOTDJylERcGiRaLXGx4uTr+1a0Pu3ML93QiIS4pj+anlBIQGcPnxZQo7FmZm/ZkvJkv1Vep7VbdZ3FuL/66Lb0zCN2+K5Ltkifh3p55PWO7UC1QdABod1N9/AXeLwzBCxcXRmg3l6vznOplhnvA6Qm+FsvH8Rvzr+mNhasG8RvOo4FyB4rmKZ8n9M5UFC4Ru+f37ULGi2PVt2zZTb/lSArR/pkiASiTZlfQkYUegnKIohwAvYO+rD1BVdSF/J2ZPT8/MHb/OCqZMge+/F4NXNWuK3ZgmadtrNCTq3/KRD2MfMvD3gXg4e7C+9XpalmqZKZKSbzqNvu7rN24Ipc6lK5JQS/1C8WH72dX/J9zyO1BuZmu+fLiORBMwV0zwGzMXGvYGRckS84RXSZGV9A/x58/wP3G0dGRA1QEUdCho/AIbFy6IlTkrKyEoU0WImeDnl2WtFX1KgEokxkJ6kvBEYCngDoQCa/QaUXZAVWH/fjFgZWMjSs+tWomys4eHoaNLNdeeXCMgNIBb0bfY+ulW3HK4cbbvWYrnKp6pAy6psfi7eRMG+Yey/a/dKJZR2IzazDPNTbAqgM3jQ+DmzZBaX+C1JImgWgXxq9z6XyXezDZPeJVLjy7RbG0zLjy8gHsOd2bWn0l3j+7Ymmd/Y403oqoQFCT0nHfsENPOvXtDnz7Qt6+ho5NIPgjknvA/iY+HNWtEz+v06ZdvSkbG0TtH8Q/xZ+P5jZiamNKxXEfmN56PmcYsS+7/ak8YXu58Vs7tyvjxsGhnCNpOH4MmARQoRj6mhVjQaM8NTPr0hR9/zJJY30ZUfBRXHl/B08WT+OR4mq9tTpcKXWhTpg2mJkasc6PTid/zgAA4cUK0VPr3FwNXeY14iEwiySZk6p7we0lyMowbB/PmiT5YuXJCXrJDB0NHlmZ+Pv0zHTd3JIdFDkb6jmSw12Dy2eVL0zUyOtn8ulNqt4ql2LMoH203XEVnE4GZ93a0miRQwEQH3fbd5aN7JTBZvMTgr/vt6NvMODSDBccXkMsqF1cHXcXS1JKdHXcaNK4Mk5QEZmaivDxrlhDcXrgQOnY0CulUieR95MM+CUdEgIuL+Hv16sJebcgQofhjBCtGICZzfz7zM3lt8tK4eGOexj9l6cml9PDogZ1F2lSL4O2n2PSUeiMiYNIkmL/9KMlV/aHURmyT82OtHcB9i69R1ETMdAq1wnsQXfJTgr+oneZ76IuLDy8y4c8JrD6zGlVVaVe2HSN9R1LRuaLBYtILN26IOYY1a4RhSK5cEBn5RjMFiUSSMeRJ+G1otaL/NWMGhIYKoY08eSAwUCj+GAlR8VHMPzafWYdncTfmLu3Ltqdx8cbksMzBEO8h6b5ueiabX0dkpOj5bryyHK3TUeh6AnsTW/rdK8qAFZcY3TCBwCLjSVBOY6GW53y+UihP49Mdd3pRVZUkXRLmGnMuPbrEhr820L9Kf4Z6D8XdwXi0vV/L4cOi5Lxhg0i27dpB3N+9+reYKUgkkqzjw0nCMTFi/2XWLLh6VQjL//CDULUCo0rAMw7N4Jv93xCTGEPdwnVZ3nw5dQr/d10nPaRlsvl1PHoEk6cmMXPznyS2bgQVE0DRMeCGE+PXRGJv/YDlnh24mKcgFjpHLCj14rlZtV4EYld684XNTAmeQt3CdRlfezyNijfi1tBb5LQyDonRt3LjBnh7Q44cMGKEENfIn9/QUUkkkld4/5OwVgsaDdy9KwwUfHzETkzLlmBqPD9+2N0wiuQsgr2FPXms89CsRDNG+I7Qe6k0NZPNr+PpU5g8PYaA/YtJqBSAU5NiPDBLRIcOjQ5cniRhP20OdO1KjktRxGw6A1m4XpRCyq701JCpXH1ylaI5i1IytzDVMFFMjDcBx8bCsmVw7ZrYWy9YUKi51a6dZjMFiUSSdbyfPWFVheBgUXJWFPjlF/H1ixehROa/0esLVVXZfXU3/iH+BF4PZNon0xjmMyxT75nWnnBMDEycdZ+AP2cTX34uWD3BJ96NTy8kM8YzikRtIuYmZgR22oOPe/V/3Ser1ov+Sbet3Vh2chlVXasyyncUzUs2z5Rd6Szj7t2XZgqPHwvJ1KAgMYAlkUgMQlp6wu9XEk5MhPXrRfI9flw4F/XtK8rORjJoBSL5rjy9kqkhUzlz/wwudi4M9hpM78q9yWGZI9Pvn5oEGR8Po2eHsjgwiOf5doL7HzR+4sZXW+/gfUuFVq0I/aYrQVGn8CvoZzAZx5tRNwkIDaCPZx9K5SnFufvneBT3iBoFahi/GcDGjWKSPClJyEqmmCkY+88lkRg5H+5g1rRp8OWXwrN3/nzhI2tt/e7nZRMStYmYa8xRFIVlJ5ehU3UsbbaUDuU6YK4xz7I43qZclJQE3y04yrTTY0jMGwLeSZgrGlb8otLuxiP4fIAo+xcqhA/gQ8Msi/ufnLx3Ev8Qf9adXYeiKFRwrkCpPKUok7eMQeLRC6oKe/cKMwUfH5Fwe/USr3fRjBluSCQSw2DcJ+GzZ8XqRdOmQkYyMlKID9SrZ1SrF7ee3mLW4VksO7WMk71P4mrvysPYh+SyypVtTmtarcoXS3Yy+8QU4p2DME02RavRoioqGkXDDzaN+aLHMnBwMGicqqrSYl0Ltl7ciq25Lb0r92aw12DccrgZNK4MkZj4Ulzj9Gkxz7Bxo6Gjkkgkb+D9PgnrdPD776LkvHevEBkoKyzQcHKCBllrPp8RTt07xdTQqaw9uxZVVWlbpi1JuiQAclvnNmhsobdCCboRRE13P+6dqMpn+72Iz3kcB1t7xh6wodKN5zTtaEKiqQnmGnP82o42WAJO1iWz5+oe6hetL069ThXwzu9NH88+OFga9kNBhlm0SJgp3L0rfs+XLoX27Q0dlUQi0RPGl4SbNIHffgNXV2G707OnEB8wMiJjIqm8sDKWppYMqDKAwd6DKehQ0NBhASIB115Rm/jkBBStBbqlgdQo50a3yNN8dioa83oN+XN0F4o+jORW3DHcrD2JfFhAuEtnIbFJsSwNW8q00Glcj7pOSPcQfNx8GFtrbNYGom+uXRMiMpaWYvItJfl+8ons90ok7xnGV47euFE0Jlu1MqoJ0CRtEuvOreN4xHGm158OwObzm/Er6IejlaOBo3tJZEwkDZd+yolHQaAAOhNaOo5jnXddTBcvgKFD2ZKQQ6+qWmklJjGGaSHTmH1kNo/iHuGd35vR1UbTtERTTBTjaUP8h9BQMdeweTP89BN07y76wDLxSiRGxftdjm7VytARpInohGh+Ov4TMw7P4Hb0bcrkKcPzxOfYmNvQolQLQ4f3ggfPHzBg09dsuLIcnZKAiQqKKiwER7Tyw9TNE7zF75T/pH16UdVKKwnJCViYWmCimDD36Fx83XwZVW0U1dyqZZveeZpRVZF0p02DkBBR0h81CurXF9831p9LIpGkCuNLwkZE4LVAWq5vSXRCNLUK1mJB4wXUL1o/W5zWUnq+lV0qU8r8E7763oLtDmtpf8GU70ISiCxRgD9aV8Kv6eD/rBdlVFUrrYTdDcM/xJ+we2Gc7XsWazNrLg28ZNz93hQRGRDmIVFRQs2tWzcx/SyRSD4IZBLWM6cjT/M88Tk+bj545POgWYlmDPIahKdLqioTWcLEvVv4KrgtqpoEKGiWH0ATUYNfy37GR1bnMV8+nGING1L9DRPm6VXVSguqqrLv+j4mB09mz7U92Jnb0btyb+KT47ExtzHeBJwirrFqFZw6JU6+27ZBvnwvk7JEIvlgkElYD6iqyt5re5kaOpXdV3fzkftHHOh6gJxWOVnRYoWhw3tBojaRETt+5McTY1FNhI2golMpXGEuI+YUo069WamS8hxZr8Rre8L6lJ3cdXUXDX5ugLOtM5NqT6K3Z2/jTbwg1ummTYPVq1+Ka0RHiyQsNZ0lkg8WmYQzyI5LO/hq31ecjjyNs60zEz6eQG/P3oYO67X8cmYjs08OxSXKgof2oFXARNVg6l6M5Wf+olejj1N1ndf5BWdUdjI2KZZlJ5ehoNC3Sl/qFq7LqharaF26NRamxmOu8VquXxce1dbWYpp/yBApriGRSACZhNPF0/inmGnMsDaz5u6zuyTrklnSdAkdynXINgkj9FYo2y9t53b0bbxcvXAK78/ELxuzVpcDrxt2fO/nza+eZpiYVCLWqhRxaeznvk1VKy08in3E3KNzmX1kNg9jH9K4eGP6VumLxkTDZ+U/y/D1DUJiIqxdC1euwPffQ6FCsHKl2GE3wnU6iUSSeRjfipIBSVG2WnhiIeNqjWOg10CSdcloFE22ms5de3YtHTd1RKtqQYW6lxwJXPOAkqU1FCq6ib+Km6J7pf/o6mBF8JjUnYT1xZKwJQz8fSCxSbE0Lt6YUb6jqF6gerZ6LdPEkyewcKEYsIqIgAoV4OhRo1qlk0gkGSctK0qGH9M1AsLuhtFxU0cKzyrM9EPTaVisIR+5fwSAqYlptkoa4w6Oo/3G9mh1ol+rUcH1rjMrZz7h9Gn4/H9eWFj+W4c6q2wEQQyu3Xp6C4ASuUrQpnQbzvY9y/b226nhbsSmCtu2CY/qMWOgVCkhKBMWJhOwRCJ5K7IcnQoG7RzEyXsnGVh1IIO9BuPu4G7okF6gU3X8dvk3KjpXJL99fsredKBrGKwtC4kmCiYmFnT5aTF+RYUMZmb0c9+FqqoE3QhiSsgUdl7ZyYAqA5jdcDbVClSjWoFqmXbfTOfoUTHRXKkSeHgITedhw6Cifj2eJRLJ+4ssR79CojaRNWfWMOfoHHa034GTrRMXH17EydYpW03nHrx5kFmHZ3E84hg3nt7kC4tGEL2DGQE6+ibP5smAErg2DKNhKcPZCAJsubCFCX9M4GjEUfLa5GWw12D6evbNViphaUKngx07xKTzwYPQrBls2WLoqCQSSTbi/VbMyiSi4qNYcGwBs47MIuJZBGXzluXOszs42TpRIrf+SrX6MLMfsnMIsw7PREWoWn0XBLVCo/FLVPnsMxMGjhtMwYIA9fUWd1pIsWQE+PXSrzyJf8L8RvPpUrELlqaWBolJL6xbJ8wULl2CAgWEq1GPHoaOSiKRGDEyCSMSsPsMd6IToqlTuA6Lmy6mXpF6eu9Pbgm786/92jtRcXyx6QzAOxPx47jH5LTKCcCBE5tRVf7e84WftV35o9oSjk9VqFRJ3OeztVlXbk4hKj6KeUfnMfPwTLZ+uhWv/F5M/WQqtua2aEyMVIjiwQPIkQPMzeHmTbCzE7aCrVunaqdaIpFI3sYHO5h1LOIY00KmAeBg6cAPtX4grHcYezq9tMTTN/67Lr5Rc/lNXHx4kc+3fU4+f2dOnfgNgIG2ozFP1oBWg6pa0XdgL/YGvkzAX2w6w52oOFReJvotYXf0/vOkcDv6NiN2j8Btuhtf7vuSis4VX6xq5bDMYZwJ+OJF6NNHnHjXrxdfGzZM9IE//VQmYIlEohc+qHeSlCGmaaHTCLoRhIOlA909uuNo5cggr0GZfv/UaC5vCbvDt79v4nrc7yimt4nhLBZa6HFCRXNjA62eNmTTpn7kruhBzS5BDGnhR3X3lz3ftyX6zDgNJ2mT8FzoycPYh7Qr245RvqOo4FxB7/fJElQV/vwTpk6F7dvF6bdLF/DyEt+XiVcikeiZD+Zd5cTdE3Tc1JHzD8/jZu/GtE+m8Xmlz7G3sM+yGN6lubwl7A5DNq0nXDMGVZMIKnQLg6/P5OeE+2iqr+xGsoXQfxg2zAcbm/8OXGWFuUJweDA/n/mZOQ3nYKYxY3HTxZTJWybb+CFniH79hL7zN99A//6QN6+hI5JIJO8x73USfhT7iMjnkZTOUxo3ezccrRz5ueXPtCndBjNN1u9vvklzeUidgiwNW8pXm38iVlcUVZMMCpjoIMysESUjNqC9aUmPHiIBOzu/+R6ZZa6gU3X8eulXJgdPJvhWMDmtcjLEewjFcxWnUfFGGbq2wYiJgSVLYOlSOHAA7O3hl19ECdra2tDRSSSSD4D3MglfeXyF6aHTWXpyKRWdKxLSI4Q8NnkI7h5s0LhSysHf/r6JW7HHyGdVhirFoxkS2JmI+AeUfwA22lIccjFFVZPRaS04efgrrFyiObHNkrJl332PzDBXuBF1g0arG/HXg79wz+HOrPqz6O7RHRtzm3Rf06BERMDs2TB/vrAQrFYNIiNFEi5Z0tDRSSSSD4j3KgkfizjGxD8nsvn8Zsw0ZnQs15FhPsMMHda/cModzhV1NAnmCUQl6zj/F9S+BkuPW/AkRwOmODfC9GhbkmxPoImqQi5fE4p6nKVs2aw1V3iW8IzzD89T1bUq+e3zU9ixMF9W/5K2ZdoapIqgN8LDhXmCVivENYYPB29vQ0clkUg+UIw+CWt1WnSqDjONGYdvH2b/9f18Uf0LBlQdQD67fIYO71+cf3Cen078RKI2EZ2qQ1Gh3zlr5nh+y60Rveg42oITm6zR2MSTs0Y+bH1uYW2RzMh65dJ0n4yYK0TGRDL7yGzmHp2LhcaC8KHhmGvM2d5+e7quZ3BUFQIDhZXgkCGi1OzvD40bQ5Eiho5OIpF84BitYtbzxOcsO7mM6YemM8J3BH08+xCfHI9Wp812ZdLg8GCm7B/Hths7KRBnzgM7jRC0UEzZ1mI3+9d8RECAeGzjz6K55RZGZFxMlu743oy6yeTgySw9uZSE5ARalGrBKN9ReOX3yvR7ZwpJSUJcY+pUOHVKJN9Ll8Aie7hcSSSS95f3WjHrXsw95hyZw7xj83gc9xgvVy8KOxYGyDZqTKG3Ql+sQK068hMhD8PIFQvfHYH+9jW4Mn44+x6f5PlZPzrV9eHePfjsM5gwAQoUsAdqZlmsWp0WjYmGq0+usjhsMZ3Ld2aE7wi9qoRlOfv2QefOcOcOlC4NixZBx44yAUskkmyH0SXhdhva8cfNP2hesjnDfYbj6+abrZx3Dtw4QIOfG5CoTUSDCTmfJjH7kCndynyGzdxRULo0Z/fDL0MbcOoU+PoK6WGvLDxwqqpK4PVApgRPoUyeMkyvP51aBWsRPiQcJ1unrAtEn9y6BbGxUKKE8O8tVQp++gnq1QOTD1aTRiKRZHPSnYQVRfkR+F1V1SxtFgZ8EoC9hT3FchXLytu+k6fxT5l/ZC7jg34gXo1HBVBgQM56DNi0HJycuHoVRrQQSbdgQVEtbdMGsuozhFanZeP5jUwJnsLxu8dxtnWmcfHGACiKYpwJOCxMmCmsXQsNGgiRjUKFYM8eQ0cmkUgk7yRdSVhRlBqAc1YnYIDKLpXT/Vx9mCe8yp3oO8w8OIX5xxfwjASq3IbTLgrJGhPMNeZ83Pk7oq2cGDcKZs4U9rITJsDQoWCZxdXzUXtGEXAogGI5i7Gw8UI6VeiUbUr4aWbfPhg/XvxpZyeGrgZlvuqZRCKR6JM0J2FFUcyAn4DfFEVppqrqVv2HpX8yYp7wz2uk7Pi6WXvyfYOWrNvbivUxh2l3DkbGV8Kj9/8IrZCToJsHqVHAj3O7fGj+lfAB6NpV5I18WTS0nWKo0LBYQyo4V6Bn5Z74uvnSvGRz49RzTkgQ/r2mphAaKvSdp0yBXr2EyYJEIpEYGek5CXcG/gKmAAMVRSmgqups/YalfzKqqbwl7A6DN60jXDMGTJN4mriGIZuSmeTaiQln81Loi++gsjil+wCJ16sxsDmcPCm0IH777cW3M52IZxFMD53OguMLeJb4DEVRqOBcgZK5S1IytxGKUTx5IoQ1Zs2CGTOgXTtRShg5Uug7SyQSiZGSnokVD2Chqqr3gFVArVcfoChKL0VRjimKcuzBgwcZjVEvZERTWafq+OLXRTzSjQUlCRSAJJ6qJ5kdVYpCy7e9yLA3bgiXOz8/ePxYtCr/+CPrEvCwXcMoNLMQAYcCaFS8EWG9wxhTfUzW3FzfXL8OgweDmxt8+SWULy9WjUDISsoELJFIjJz0nISvAIX//rsncPPVB6iquhBYCGJPON3R/YOM9nPTq6msU3VUmViYC9qb5IuHBA1oFQUUMyx15V4k8efPYdIkoQOh0QiN5xEjwCpjks2pIuxuGBWdK6IoCnbmdnzu8TnDfYe/WN0ySlRVCGpcugQdOghlq/LlDR2VRCKR6JX0JOHFwBJFUT4FzIDW+g3pv+ijn5saTeWUnm94bAg5NUlMazaDFpXy0+m+M5+dSOSqQwe2lM5LrOk5LHXlsNCVIl8OK9asgVGj4PZtaN9etCnz59fjC/AaVFVl55WdTA6ezIGbB/j9s9+pX7Q+Y2uNzdwbZxY6najZL1wIa9aAjQ0sXixOwa6ZL1YikUgkhiDNSVhV1WdAm0yI5Y3owyP3XZrKW8LuMGDTAiI041FNdTxVYMJSCxRlCkMm7mfL+Uf8tPksZklaciSXAUB56MiTXZXpcBI8PETuqF5djz/4a0jWJbP+3HqmBE/hVOQpXO1cCfgkgGpu1TL3xplFfDysWiXWjC5cEEn30iXxgkpNZ4lE8p5jFGId+vLIfZOm8oOY+3y5phkRVsdf7PcqOrhvkfB3ov+Y5pXyg6Lgv+sit+5oSThchgfH8pE7t8JPP0G3bqIMnVmoqoqiKOhUHYN+G05sgjm5EofgllCfQpZlsbOwy7ybZxaRkVChgvjTwwN+/lksTpsZsUGERCKRpAGjkBJ6U982Ix65qqpy//l9AGyTFGITT9DkkimmOhNQTUAxJ9Hc+1+JvlFZV9pafUz0yro8DnNhyBCFS5fg888zLwE/iXvCuIPj8PzJk0RtIr+dfkCOmInkiZ2DrbYOd58m88WmM2wJu5M5Aeiba9dEyQDAyQk6dRIGC8ePi96vTMASieQDwihOwvr0yP0z/E/mB88k7NJB4uKfcWncU6wc81AreRYhrq7kTr5KvMmZFz3flES/fz8MHAjnzsEnn4hNmVKl9PYj/ofb0bdfrBk9T3pOw2INeRT7CP9dF9Em5eGfIltpLc0bhCNHhJnCxo2i39u0qfjT39/QkUkkEonBMIokrA+P3PjkeP63ZQhTzi5AVUBRYdjDAugePYS8+WjWowWHNp3BIqkUFjqRXa3MNHQtX5q2beGXX4TU5ObN0KzZm6Um9aHKdfb+WSotqIRO1fFp2U8ZVW0U5Z3EZHBE1InXPietpfks4+RJoWT1xx9CUGPUKPFpxiZ7OV1JJBKJITCKJAwZ88gF2L1tOpPPLXjxbxNFQ67OfTDPm+/F9eFlone2saHwvUoMbm2PqsLYsUIb4m0rRxmZ4j58+zCXH1+mY/mOlMlThu9qfsdn5T+joEPBfz0uvatWWUp8vBDYyJdPJNvbt2H6dOjRQ0hMSiQSiQQwYj/hd3Hr8Q1mrB6EEzaMGrAGXWICCya1ZrjJXhJ1SZhrzAnsHIiPm8+/nqeqwgNg6FDRvmzVSgzuuru/+57VJu17bYJ0dbAieMzH//m6qqrsurqLSX9O4sDNA7jncOfKoCuYmrz5s9GriR7EiX1iy3KGL0c/egTz5sHs2cIeavNm8XWdTjoZSSSSD4b32k/4Xaw6soiA3f/jdJIYVOp3W5x0Tcwt6Pvtdir+7fXrV9DvPwn4yhVROf39d9Hv3bsXatdO/b3TMsUdHB5M/9/6cyryFPnt8xPwSQA9K/d8awIG/ZTm9c61axAQAEuWQFwcNGwolK5SkAlYIpFIXst7lYR7B9Ri4bMgUMUPts59JC2/nfivx/i4+fwn+cbGwsSJQmTDwkKcfAcOTPuSb+FwAAAPzUlEQVSg7rtKxXFJccQkxpDHJg/WZtYkahNZ2mwpHcp1wFyTegnGjJbm9Yaqiub48uXCu7djRxg2DMqUMXRkEolEYhQY9RFFq9Oy8bdpXL0QAoDGxg5FBRRQNRouFnN86+6QqoqKaenSMG4ctG0rjHmGDUvfpszIeiWwMvv3/azMNPT72JkJf0zAfYY7I/aMAMAjnwfn+p2ja8WuaUrABkenE/X6jz4SxsggbARv3BAKVzIBSyQSSaoxqpNw6N+lZB9Xby7+sYmp5xdzxSqOL4K9GD8ulI4NxrBs5V4StYmYa8zxK+j3xmtduiRKz7t2QblycOCAyCsZ4dVSce4cz3F120efwJ+JSYyhQdEG9PDo8eLxyptGrLMjrypbFSgA2r/70o6Oho1NIpFIjBSjScKht0KpvaI28clxqH+fdqs8M2OD/Wc0HzADRVHwLeBLYOfAN/Z8QRgtjB8vcomlpdj37d9fWNTqg3+WiofuHMqsIz/Rrkw7RlcbTQXnCvq5iSGoUweCg4Wy1erVQtlKXy+aRCKRfKAYzbto0I0gErWJpMxy97TzY8GYnSgWFi8eI3Z044iIKs8OhzhG1rvzIiGqqqieDhkC4eHQuTNMngzOzvqL8eido0wOnkz/Kv2pVagWo6uPZqDXQON0M7p2TXj4jh0r9rK++EL8WavWm5ekJRKJRJImjCYJ+xX0w1xjLkrNpuZ0azPhPwn4TTu65exdGThQTD2XKyd0I/RltKCqKnuv7WVS8CT2Xd+Hg6UDTUs0BcDZVo8ZPqs4elSoWG3cKPrpDRqIxNuokaEjk0gkkvcOo0nCPm4+by01v85pKTZOZdDIBO7/KQatAgLE1LM+q6hN1jTh18u/4mLngn9df3pV7oW9hb3+bpBVPH0qpCQPHhTKViNHihdL2ghKJBJJpmE0SRhev16Uwqu7uHFX8/B4bxmSo2z49FPRA3ZxyXgM8cnxrD+3ng7lOmBqYkrLUi1pUbIFHct3xMLU4t0XyE4kJMCJE+DjA/b2kDu3eKE+/1z8WyKRSCSZilEl4beRsqObHG3Jk8DSxF7Kh2nOGEr3CGPNIo8MXz86IZr5x+Yz/dB07sXcw9HSkSYlmtDdo7seos9injwR/d5Zs8QJ+PZtyJlTlKAlEolEkmW8N0l46Mcl6PflMx4eLAqqgsNHF8jre5Pxbctm6LpxSXGMOziOuUfn8jThKXUK12Fli5XULpQGKa3sQkSEUCRZtEiMidetKwwV5IqRRCKRGIT3IgkHB8PYPq48OAuOpR5gXfMM7u4wsl7ZdCtLxSTGYGtui4WpBdsubaNO4TqMqT4GT5dUyYFmLxITwdxcnHp//BHatYMRI6CCEa9MSSQSyXuAUSfhR49g9Ggh1OTmJlaQmjXLA/zXLCG1nI48zeTgyey5uoerg65iZ2HHkc+PYGWWjVyKUoOqCiUSf3/IlQvWrxeC2BERovcrkUgkEoNjlLKVqgpLl0KJErBsmRjk/esv4fObvuup/HHzDxqtbkSF+RXYdnEbXSp0IVmXDGBcCTgxEVasEKfcBg2EDqe398vvywQskUgk2QajOwmfOwd9+4pd32rVhHNeuXIZu+aZ+2f4aNlH5LbOzbha4+hXpR+OVkbaJ50yBb75BsqWFZ9Q2rcXpWiJRCKRZDuMKgmvXAndu4vtmUWLoFu39LnkJWmTWHduHeFPw/myxpeUdyrP+tbraVS8EdZm1voPPDO5fRtmzhRDVp98ItaLPD2hXj2pbCWRSCTZHKMqR1erJhLvxYvQo0faE3BsUixzjsyh2OxidNrciU3nN70oObcp08a4EvCZM9ClCxQqBNOnw7Fj4uvOzlC/vkzAEolEYgQY1Um4cGFYuDB9z915ZSedN3fmQewDfN18mdNwDg2LNcREMarPIYKePUUpwMYG+vWDoUOhYEFDRyWRSCSSNGJUSTitRDyLIDYplqI5i1IsZzGqulZldLXR1HCvYejQ0kZyMmzaJCbPLCzA11ck3b59hciGRCKRSIwSRVXVdz8qA3h6eqrHUkqlWcTlR5eZEjyFFadXUK9IPba135al99cbMTFi/2r6dLh5U1gItm9v6KgkEolE8hYURTmuqmqqRCXeq5PwyXsnmfDHBDb8tQFzjTk9PHowwneEocNKOwkJ8P33YvT7yRPRDJ85E5o0MXRkEolEItEjRp+EU07yiqKw9cJWdl3dxehqoxniPQQnWycDR5dGHj8W5WVzcyG04ecnlqB9Xm9aIZFIJBLjxmjL0TpVx7aL25j05yRG+o6kVelWRCdEo6oqOSxz6P1+mYaqCt1Nf38IChJlZwcHcRq2MDJXJolEIpG83+XoJG0Sq8+sZnLwZM4/PE8hh0IvJpyNysdXq4WtW0XyPXRInIAHD375fZmAJRKJ5L3H6JJwvVX12H9jP+WdyrO65WralGmDqYnR/Rhw/jy0aiX2fGfPFgvQNjaGjkoikUgkWYjRlaO3X9yOxkRDg6INUIxJkOLRI+Fg9PixmHYG2L8fatQAUyP8ECGRSCSS1/Jel6OblDCyCeHr1yEgAJYsgdhYseur0wm5r1q1DB2dRCKRSAyI0SVho2LlSujaFTQa6NBBePiWLWvoqCQSiUSSTUi3ZqOiKE6KooTpMxijR1Vh5044ckT8288Phg8Xp+Fly2QClkgkEsm/yIhw8lTAiIx2M5FXPXxTer5ubsJa0NXVsPFJJBKJJFuSriSsKMrHwHPgnn7DMUKWLIEiRYSjkaqKE+/y5YaOSiKRSCRGQJp7woqimAPfAC2ALXqPyBiIiIA8ecDMDB4+hKJFhb2TtBCUSCQSSRpIz0l4DPCjqqpRb3qAoii9FEU5pijKsQcPHqQ/uuzGX39B9+7CwWjdOvG1ESPEqlGDBjIBSyQSiSRNpCcJ1wH6K4oSBFRUFGXRqw9QVXWhqqqeqqp65smTJ6MxGhZVhYMHhXlCmTKwdi306iVMFUCsGkkkEolEkg7SXI5WVfWjlL8rihKkqurn+g0pG9KvH0RGwtix4u+5cxs6IolEIpG8B2RoT1hVVT89xZF9iIsTk86LFkFgINjbw8aNYtLZ2trQ0UkkEonkPULWUlN49Ah++AHc3aFPH/G1iAjxZ4kSMgFLJBKJRO9IxSyAO3egeHEhK9mwofDwrVlTDlr9v737C5GrPOM4/v0ZrVjbptrapQ2NRKgX0kRMq3Vp1VVMDGmgQYqKLQSkIKEGexFQSUWEVsiFQgmNxLYXbaViQ+lFtZIUJCj2jxuN0YoWbUms9W9pqc1VwDy9OGdxnJyzebvunvc9Z38fWDK785zN85xn533mzJyZMTOzBbV4h/CBAzA9XZ1ktWwZ3H47bNjgd7UyM7POLK6HoyNg715YswZWr4Zt26qjX4Bbb/UANjOzTi2eITw9DRdcAFddVb3ed/t2ePllP9drZmbZDPvh6CNHqs/vXb68elnRsWPV20xefz2cemru7MzMbJEb5hB+803YsQN27oTJSXj4YVixAg4e9MlWZmZWjGEN4Zdegrvvrj5E4ehR2LixOtN5hgewmZkVZBhDOKIasA8+WA3gTZuqz/E999zcmZmZmbXq74lZx47BQw/BpZfC7t3Vz7ZsgcOHYdcuD2AzMyte/4bw0aPV0e7KldWHKhw6VB0JAyxdChMTObMzMzNL1r+Ho9etqz46cNUquP9+uOaa6nN9zczMeqZ/Q3jrVrjlFli71idamZlZr/VvCK9fnzsDMzOzedG/54TNzMwGwkPYzMwsEw9hMzOzTDyEzczMMvEQNjMzy8RD2MzMLBMPYTMzs0w8hM3MzDLxEDYzM8vEQ9jMzCwTD2EzM7NMPITNzMwy8RA2MzPLRBGxsP+B9DZweB5/5SeBf87j78vJtZRpKLUMpQ5wLSUaSh0w/7WcHRFnpQQu+BCeb5L2R8QXc+cxH1xLmYZSy1DqANdSoqHUAXlr8cPRZmZmmXgIm5mZZdLHIXxf7gTmkWsp01BqGUod4FpKNJQ6IGMtvXtO2MzMbCj6eCRsZmY2CMUOYUk/kfQHSd/9IDG5SVoq6RFJeyX9WtKHGmJOlvSKpH3118ocuc4mNUdJd0qalvTDrnNMJWnzSB3PSNrVENOHnkxIery+fIqk30h6QtINs2yTFNe1sVqW1/v8UUn3SVLLNsskvTrSo6SXhCy0sVqScyxtPRur486RGl6UdFvLNsX1pGkNTt3XXfSkyCEs6WpgSURMAudI+txcYgrxDeCeiFgLvAGsa4hZBTwQEVP113OdZpjmhDlK+gLwFeAi4C1JV3adZIqIuHemDuBx4EcNYUX3RNIZwE+B0+sfbQGeiogvA1+X9NGWTVPjOtNQy43A5oi4Avgs0HYH6EvA90d69PbCZzu7hlqScixtPRuvIyLuGLnN/Bn4WcumxfWE49fg60jY1131pMghDEwBv6wv76Va2OcSk11E7IyI39XfngW81RB2MbBB0pP1Pa+Tu8swWUqOlwG/iupEgz3AJZ1m+H+StAyYiIj9DVeX3pN3gWuBd+rvp3jv9vAY0Paax9S4Lr2vlojYFhEv1Nd9gvY3UbgY+JakpyXdtfBpJhnvS2qOU5S1no3XAYCkC4FXI+IfLdsV15OGNfibpO3rqcS4D6TUIXw6MNPkfwETc4wphqRJ4IyI+GPD1dPAlRFxEXAKsL7T5NKk5NirngDfBu5tua7onkTEOxHxn5Efpe774nrUUAsAkq4Fno+I11o2fYRqobwQmJS0auGyTNNQS2qORfWlrSfAzcCOWTYtriczZtZg4O8UdFspdQgfAU6rL3+E5jxTYoog6UyqP9y25+CejYjX68v7gRIfWk/JsU89OQm4HNjXEtKHnoxK3fe96JGkc4CtwHdmCft9RPw3It4FDlBmj1JzLL4vkj4OfCoi/jpLWJE9GVuDi7qtFNfo2lO8d+h/PnBojjHZ1Sdi7QZui4i299D+uaTzJS0BNgIHO0swXUqOvehJ7RLgT9H+Gr0+9GRU6r4vvkf185EPADe0HI3N2CPp05I+DKyleq6yNKk5Ft8X4GvAb08QU1xPGtbgsm4rEVHcF/AxqkXvHuCFegd87wQxS3Pn3VLLZuDfVEdc+4A7Gmr5PPAs8BzVSQ3Z826o4305AmcCPx6LOQl4AvgB8BdgRe68Z6nnLuDq+vJ5fexJnee++t+zgefrfT8NLAGuAG4aiz8uLncNDbVsB14fuc1c1lLL5cCLdZ9u6jrfxFqOy7Hl763I9WymjvryL4DVI9/3oicNa/Cm8X2dsyfFvllHfW94DfBYRLwx1xjrlqTTgK8CT0fE33Lns5hI+gzVPfc9McsRZGqcdcvrWXdS93UXPSl2CJuZmQ1dqc8Jm5mZDZ6HsJmZWSYewmZmZpl4CJuZmWXiIWxmZpaJh7CZmVkm/wPy1LuZ2QKh4QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "prstd, iv_l, iv_u = wls_prediction_std(res_wls)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(x, y, 'o', label=\"Data\")\n",
    "ax.plot(x, y_true, 'b-', label=\"True\")\n",
    "# OLS\n",
    "ax.plot(x, res_ols.fittedvalues, 'r--')\n",
    "ax.plot(x, iv_u_ols, 'r--', label=\"OLS\")\n",
    "ax.plot(x, iv_l_ols, 'r--')\n",
    "# WLS\n",
    "ax.plot(x, res_wls.fittedvalues, 'g--.')\n",
    "ax.plot(x, iv_u, 'g--', label=\"WLS\")\n",
    "ax.plot(x, iv_l, 'g--')\n",
    "ax.legend(loc=\"best\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feasible Weighted Least Squares (2-stage FWLS)\n",
    "\n",
    "Like `w`, `w_est` is proportional to the standard deviation, and so must be squared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            WLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.931\n",
      "Model:                            WLS   Adj. R-squared:                  0.929\n",
      "Method:                 Least Squares   F-statistic:                     646.7\n",
      "Date:                Tue, 24 Dec 2019   Prob (F-statistic):           1.66e-29\n",
      "Time:                        15:06:39   Log-Likelihood:                -50.716\n",
      "No. Observations:                  50   AIC:                             105.4\n",
      "Df Residuals:                      48   BIC:                             109.3\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          5.2363      0.135     38.720      0.000       4.964       5.508\n",
      "x1             0.4492      0.018     25.431      0.000       0.414       0.485\n",
      "==============================================================================\n",
      "Omnibus:                        0.247   Durbin-Watson:                   2.343\n",
      "Prob(Omnibus):                  0.884   Jarque-Bera (JB):                0.179\n",
      "Skew:                          -0.136   Prob(JB):                        0.915\n",
      "Kurtosis:                       2.893   Cond. No.                         14.3\n",
      "==============================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "resid1 = res_ols.resid[w==1.]\n",
    "var1 = resid1.var(ddof=int(res_ols.df_model)+1)\n",
    "resid2 = res_ols.resid[w!=1.]\n",
    "var2 = resid2.var(ddof=int(res_ols.df_model)+1)\n",
    "w_est = w.copy()\n",
    "w_est[w!=1.] = np.sqrt(var2) / np.sqrt(var1)\n",
    "res_fwls = sm.WLS(y, X, 1./((w_est ** 2))).fit()\n",
    "print(res_fwls.summary())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
